|
LR-splines
0.5
|
00001 #include "LRSpline/LRSplineSurface.h" 00002 #include "LRSpline/Basisfunction.h" 00003 #include "LRSpline/Meshline.h" 00004 #include "LRSpline/Element.h" 00005 #include "LRSpline/Profiler.h" 00006 00007 #include <set> 00008 #include <algorithm> 00009 #include <cstdio> 00010 #include <cstdlib> 00011 #include <fstream> 00012 #include <iostream> 00013 #include <iomanip> 00014 #include <boost/rational.hpp> 00015 #include <cfloat> 00016 #include <cmath> 00017 00018 typedef unsigned int uint; 00019 00024 typedef std::pair<double,int> IndexDouble; 00025 00026 namespace LR { 00027 00028 #define DOUBLE_TOL 1e-14 00029 00030 00031 /************************************************************************************************************************/ 00034 LRSplineSurface::LRSplineSurface() { 00035 #ifdef TIME_LRSPLINE 00036 PROFILE("Constructor"); 00037 #endif 00038 order_.resize(2); 00039 start_.resize(2); 00040 end_.resize(2); 00041 rational_ = false; 00042 dim_ = 0; 00043 order_[0] = 0; 00044 order_[1] = 0; 00045 start_[0] = 0; 00046 start_[1] = 0; 00047 end_[0] = 0; 00048 end_[1] = 0; 00049 meshline_ = std::vector<Meshline*>(0); 00050 element_ = std::vector<Element*>(0); 00051 00052 initMeta(); 00053 } 00054 00055 #ifdef HAS_GOTOOLS 00056 /************************************************************************************************************************/ 00060 LRSplineSurface::LRSplineSurface(Go::SplineSurface *surf) { 00061 #ifdef TIME_LRSPLINE 00062 PROFILE("Constructor"); 00063 #endif 00064 initMeta(); 00065 initCore(surf->numCoefs_u(), surf->numCoefs_v(), 00066 surf->order_u(), surf->order_v(), 00067 surf->basis_u().begin(), surf->basis_v().begin(), 00068 surf->ctrl_begin(), surf->dimension(), surf->rational() ); 00069 } 00070 #endif 00071 00072 00073 /************************************************************************************************************************/ 00085 LRSplineSurface::LRSplineSurface(int n1, int n2, int order_u, int order_v, double *knot_u, double *knot_v, double *coef, int dim, bool rational) { 00086 #ifdef TIME_LRSPLINE 00087 PROFILE("Constructor"); 00088 #endif 00089 initMeta(); 00090 initCore(n1,n2,order_u, order_v, knot_u, knot_v, coef, dim, rational); 00091 } 00092 00093 /************************************************************************************************************************/ 00102 LRSplineSurface::LRSplineSurface(int n1, int n2, int order_u, int order_v, double *knot_u, double *knot_v) { 00103 #ifdef TIME_LRSPLINE 00104 PROFILE("Constructor"); 00105 #endif 00106 initMeta(); 00107 00108 std::vector<double> grevU = LRSpline::getGrevillePoints(knot_u, knot_u + n1 + order_u); 00109 std::vector<double> grevV = LRSpline::getGrevillePoints(knot_v, knot_v + n2 + order_v); 00110 double coef[n1*n2*2]; 00111 int k=0; 00112 for(int j=0; j<n2; j++) 00113 for(int i=0; i<n1; i++) { 00114 coef[k++] = grevU[i]; 00115 coef[k++] = grevV[j]; 00116 } 00117 00118 initCore(n1,n2,order_u, order_v, knot_u, knot_v, coef, 2, false); 00119 } 00120 00121 /************************************************************************************************************************/ 00128 LRSplineSurface::LRSplineSurface(int n1, int n2, int order_u, int order_v) { 00129 #ifdef TIME_LRSPLINE 00130 PROFILE("Constructor"); 00131 #endif 00132 initMeta(); 00133 00134 std::vector<double> knot_u = LRSpline::getUniformKnotVector(n1, order_u); 00135 std::vector<double> knot_v = LRSpline::getUniformKnotVector(n2, order_v); 00136 std::vector<double> grev_u = LRSpline::getGrevillePoints(knot_u.begin(), knot_u.end()); 00137 std::vector<double> grev_v = LRSpline::getGrevillePoints(knot_v.begin(), knot_v.end()); 00138 double coef[n1*n2*2]; 00139 int k=0; 00140 for(int j=0; j<n2; j++) 00141 for(int i=0; i<n1; i++) { 00142 coef[k++] = grev_u[i]; 00143 coef[k++] = grev_v[j]; 00144 } 00145 00146 // generate the uniform knot vector 00147 initCore(n1,n2,order_u, order_v, knot_u.begin(), knot_v.begin(), coef, 2, false); 00148 } 00149 00150 void LRSplineSurface::initMeta() { 00151 maxTjoints_ = -1; 00152 doCloseGaps_ = true; 00153 maxAspectRatio_ = 2.0; 00154 doAspectRatioFix_ = false; 00155 refStrat_ = LR_FULLSPAN; 00156 refKnotlineMult_ = 1; 00157 symmetry_ = 1; 00158 element_red = 0.5; 00159 element_green = 0.5; 00160 element_blue = 0.5; 00161 basis_red = 1.0; 00162 basis_green = 0.83; 00163 basis_blue = 0.0; 00164 selected_basis_red = 1.0; 00165 selected_basis_green = 0.2; 00166 selected_basis_blue = 0.05; 00167 lastElementEvaluation = -1; 00168 } 00169 00170 /************************************************************************************************************************/ 00175 LRSplineSurface::~LRSplineSurface() { 00176 for(Basisfunction* b : basis_) 00177 delete b; 00178 for(uint i=0; i<meshline_.size(); i++) 00179 delete meshline_[i]; 00180 for(uint i=0; i<element_.size(); i++) 00181 delete element_[i]; 00182 } 00183 00184 00185 /************************************************************************************************************************/ 00189 LRSplineSurface* LRSplineSurface::copy() const { 00190 generateIDs(); 00191 00192 std::vector<Basisfunction*> basisVector; 00193 00194 // flat list to make it quicker to update pointers from Basisfunction to Element and back again 00195 LRSplineSurface *returnvalue = new LR::LRSplineSurface(); 00196 00197 for(Basisfunction* b : basis_) { 00198 Basisfunction *newB = b->copy(); 00199 returnvalue -> basis_.insert(newB); 00200 basisVector.push_back(newB); 00201 } 00202 00203 for(Element *e : element_) { 00204 Element* newEl = e->copy(); 00205 returnvalue -> element_.push_back(newEl); 00206 newEl->updateBasisPointers(basisVector); 00207 } 00208 00209 for(Meshline *m : meshline_) 00210 returnvalue -> meshline_.push_back(m->copy()); 00211 00212 returnvalue->rational_ = this->rational_; 00213 returnvalue->dim_ = this->dim_; 00214 returnvalue->order_[0] = this->order_[0]; 00215 returnvalue->order_[1] = this->order_[1]; 00216 returnvalue->start_[0] = this->start_[0]; 00217 returnvalue->start_[1] = this->start_[1]; 00218 returnvalue->end_[0] = this->end_[0]; 00219 returnvalue->end_[1] = this->end_[1]; 00220 returnvalue->maxTjoints_ = this->maxTjoints_; 00221 returnvalue->doCloseGaps_ = this->doCloseGaps_; 00222 returnvalue->doAspectRatioFix_ = this->doAspectRatioFix_; 00223 returnvalue->maxAspectRatio_ = this->maxAspectRatio_; 00224 00225 return returnvalue; 00226 } 00227 00228 00229 #ifdef HAS_GOTOOLS 00230 /************************************************************************************************************************/ 00239 void LRSplineSurface::point(Go::Point &pt, double u, double v, int iEl, bool u_from_right, bool v_from_right) const { 00240 std::vector<std::vector<double> > res; 00241 point(res, u, v, 0, iEl, u_from_right, v_from_right); 00242 pt = Go::Point(res[0].begin(), res[0].end()); 00243 return ; 00244 00245 } 00246 00247 /************************************************************************************************************************/ 00254 void LRSplineSurface::point(Go::Point &pt, double u, double v, int iEl) const { 00255 std::vector<std::vector<double> > res; 00256 point(res, u, v, 0, iEl); 00257 pt = Go::Point(res[0].begin(), res[0].end()); 00258 return ; 00259 00260 } 00261 00262 /************************************************************************************************************************/ 00272 void LRSplineSurface::point(std::vector<Go::Point> &pts, double u, double v, int derivs, int iEl) const { 00273 00274 std::vector<std::vector<double> > res; 00275 point(res, u, v, derivs, iEl); 00276 00277 // clear and resize output array (optimization may consider this an outside task) 00278 pts.resize((derivs+1)*(derivs+2)/2); 00279 for(uint i=0; i<pts.size(); i++) { 00280 pts[i].resize(dim_); 00281 pts[i].setValue(0.0); 00282 for(int j=0; j<dim_; j++) 00283 pts[i][j] = res[i][j]; 00284 } 00285 return; 00286 00287 } 00288 #endif 00289 00290 /************************************************************************************************************************/ 00299 void LRSplineSurface::point(std::vector<double> &pt, double u, double v, int iEl, bool u_from_right, bool v_from_right) const { 00300 std::vector<std::vector<double> > res; 00301 point(res, u, v, 0, u_from_right, v_from_right, iEl); 00302 pt = res[0]; 00303 } 00304 00305 /************************************************************************************************************************/ 00312 void LRSplineSurface::point(std::vector<double> &pt, double u, double v, int iEl) const { 00313 std::vector<std::vector<double> > res; 00314 point(res, u, v, 0, iEl); 00315 pt = res[0]; 00316 } 00317 00318 /************************************************************************************************************************/ 00328 void LRSplineSurface::point(std::vector<std::vector<double> > &pts, double u, double v, int derivs, int iEl) const { 00329 point(pts, u, v, derivs, u!=end_[0], v!=end_[1], iEl); 00330 } 00331 00332 /************************************************************************************************************************/ 00342 void LRSplineSurface::point(std::vector<std::vector<double> > &pts, double u, double v, int derivs, bool u_from_right, bool v_from_right, int iEl) const { 00343 #ifdef TIME_LRSPLINE 00344 PROFILE("Point()"); 00345 #endif 00346 std::vector<double> basis_ev; 00347 00348 // clear and resize output array (optimization may consider this an outside task) 00349 pts.resize((derivs+1)*(derivs+2)/2); 00350 for(uint i=0; i<pts.size(); i++) 00351 pts[i].resize(dim_, 0); 00352 00353 if(iEl == -1) 00354 iEl = getElementContaining(u,v); 00355 for(Basisfunction* b : element_[iEl]->support() ) { 00356 b->evaluate(basis_ev, u,v, derivs, u_from_right, v_from_right); 00357 for(uint i=0; i<pts.size(); i++) 00358 for(int j=0; j<dim_; j++) 00359 pts[i][j] += basis_ev[i]*b->cp(j); 00360 } 00361 } 00362 00363 #ifdef HAS_GOTOOLS 00364 /************************************************************************************************************************/ 00374 void LRSplineSurface::computeBasis (double param_u, double param_v, Go::BasisDerivsSf2 & result, int iEl ) const { 00375 #ifdef TIME_LRSPLINE 00376 PROFILE("computeBasis()"); 00377 #endif 00378 std::vector<double> values; 00379 HashSet_const_iterator<Basisfunction*> it, itStop, itStart; 00380 itStart = (iEl<0) ? basis_.begin() : element_[iEl]->constSupportBegin(); 00381 itStop = (iEl<0) ? basis_.end() : element_[iEl]->constSupportEnd(); 00382 int nPts= (iEl<0) ? basis_.size() : element_[iEl]->nBasisFunctions(); 00383 result.prepareDerivs(param_u, param_v, 0, -1, nPts); 00384 00385 int i=0; 00386 //element_[i]->write(std::cout); 00387 00388 for(it=itStart; it!=itStop; ++it, ++i) { 00389 (*it)->evaluate(values, param_u, param_v, 2, param_u!=end_[0], param_v!=end_[1]); 00390 00391 result.basisValues[i] = values[0]; 00392 result.basisDerivs_u[i] = values[1]; 00393 result.basisDerivs_v[i] = values[2]; 00394 result.basisDerivs_uu[i] = values[3]; 00395 result.basisDerivs_uv[i] = values[4]; 00396 result.basisDerivs_vv[i] = values[5]; 00397 } 00398 } 00399 00400 /************************************************************************************************************************/ 00410 void LRSplineSurface::computeBasis (double param_u, double param_v, Go::BasisDerivsSf & result, int iEl ) const { 00411 #ifdef TIME_LRSPLINE 00412 PROFILE("computeBasis()"); 00413 #endif 00414 std::vector<double> values; 00415 HashSet_const_iterator<Basisfunction*> it, itStop, itStart; 00416 itStart = (iEl<0) ? basis_.begin() : element_[iEl]->constSupportBegin(); 00417 itStop = (iEl<0) ? basis_.end() : element_[iEl]->constSupportEnd(); 00418 int nPts= (iEl<0) ? basis_.size() : element_[iEl]->nBasisFunctions(); 00419 result.prepareDerivs(param_u, param_v, 0, -1, nPts); 00420 00421 int i=0; 00422 for(it=itStart; it!=itStop; ++it, ++i) { 00423 (*it)->evaluate(values, param_u, param_v, 1, param_u!=end_[0], param_v!=end_[1]); 00424 00425 result.basisValues[i] = values[0]; 00426 result.basisDerivs_u[i] = values[1]; 00427 result.basisDerivs_v[i] = values[2]; 00428 } 00429 } 00430 00431 /************************************************************************************************************************/ 00441 void LRSplineSurface::computeBasis(double param_u, double param_v, Go::BasisPtsSf & result, int iEl ) const { 00442 #ifdef TIME_LRSPLINE 00443 PROFILE("computeBasis()"); 00444 #endif 00445 HashSet_const_iterator<Basisfunction*> it, itStop, itStart; 00446 itStart = (iEl<0) ? basis_.begin() : element_[iEl]->constSupportBegin(); 00447 itStop = (iEl<0) ? basis_.end() : element_[iEl]->constSupportEnd(); 00448 int nPts= (iEl<0) ? basis_.size() : element_[iEl]->nBasisFunctions(); 00449 00450 result.preparePts(param_u, param_v, 0, -1, nPts); 00451 int i=0; 00452 for(it=itStart; it!=itStop; ++it, ++i) 00453 result.basisValues[i] = (*it)->evaluate(param_u, param_v, param_u!=end_[0], param_v!=end_[1]); 00454 } 00455 #endif 00456 00457 /************************************************************************************************************************/ 00469 void LRSplineSurface::computeBasis (double param_u, 00470 double param_v, 00471 std::vector<std::vector<double> >& result, 00472 int derivs, 00473 int iEl ) const 00474 { 00475 #ifdef TIME_LRSPLINE 00476 PROFILE("computeBasis()"); 00477 #endif 00478 result.clear(); 00479 HashSet_const_iterator<Basisfunction*> it, itStop, itStart; 00480 itStart = (iEl<0) ? basis_.begin() : element_[iEl]->constSupportBegin(); 00481 itStop = (iEl<0) ? basis_.end() : element_[iEl]->constSupportEnd(); 00482 int nPts= (iEl<0) ? basis_.size() : element_[iEl]->nBasisFunctions(); 00483 00484 result.resize(nPts); 00485 int i=0; 00486 for(it=itStart; it!=itStop; ++it, ++i) 00487 (*it)->evaluate(result[i], param_u, param_v, derivs, param_u!=end_[0], param_v!=end_[1]); 00488 00489 } 00490 00491 /************************************************************************************************************************/ 00498 int LRSplineSurface::getElementContaining(double u, double v) const { 00499 if(lastElementEvaluation >= 0) 00500 if(element_[lastElementEvaluation]->umin() <= u && element_[lastElementEvaluation]->vmin() <= v) 00501 if((u < element_[lastElementEvaluation]->umax() || (u == end_[0] && u <= element_[lastElementEvaluation]->umax())) && 00502 (v < element_[lastElementEvaluation]->vmax() || (v == end_[1] && v <= element_[lastElementEvaluation]->vmax()))) 00503 return lastElementEvaluation; 00504 for(uint i=0; i<element_.size(); i++) 00505 if(element_[i]->umin() <= u && element_[i]->vmin() <= v) 00506 if((u < element_[i]->umax() || (u == end_[0] && u <= element_[i]->umax())) && 00507 (v < element_[i]->vmax() || (v == end_[1] && v <= element_[i]->vmax()))) { 00508 lastElementEvaluation = i; 00509 return i; 00510 } 00511 00512 return -1; 00513 } 00514 00515 /************************************************************************************************************************/ 00523 void LRSplineSurface::getMinspanLines(int iEl, std::vector<Meshline*>& lines) { 00524 Element *e = element_[iEl]; 00525 double umin = e->umin(); 00526 double umax = e->umax(); 00527 double vmin = e->vmin(); 00528 double vmax = e->vmax(); 00529 double min_du = DBL_MAX; 00530 double min_dv = DBL_MAX; 00531 int best_startI = order_[0]+2; 00532 int best_stopI = order_[0]+2; 00533 int best_startJ = order_[1]+2; 00534 int best_stopJ = order_[1]+2; 00535 bool only_insert_span_u_line = (vmax-vmin) >= maxAspectRatio_*(umax-umin); 00536 bool only_insert_span_v_line = (umax-umin) >= maxAspectRatio_*(vmax-vmin); 00537 // loop over all supported B-splines and choose the minimum one 00538 for(Basisfunction* b : e->support()) { 00539 double lowu = b->getParmin(0); 00540 double highu = b->getParmax(0); 00541 double lowv = b->getParmin(1); 00542 double highv = b->getParmax(1); 00543 double du = highu - lowu; 00544 double dv = highv - lowv; 00545 int startI=0; 00546 int stopI=0; 00547 int startJ=0; 00548 int stopJ=0; 00549 while((*b)[0][startI] <= e->umin()) 00550 startI++; 00551 while((*b)[0][stopI] < e->umax()) 00552 stopI++; 00553 while((*b)[1][startJ] <= e->vmin()) 00554 startJ++; 00555 while((*b)[1][stopJ] < e->vmax()) 00556 stopJ++; 00557 00558 // min_du is defined as the minimum TOTAL knot span (of an entire basis function) 00559 bool fixU = false; 00560 int delta_startI = abs(startI - (order_[0]+1)/2); 00561 int delta_stopI = abs(stopI - (order_[0]+1)/2); 00562 if( du < min_du ) 00563 fixU = true; 00564 if( du == min_du && delta_startI <= best_startI && delta_stopI <= best_stopI ) 00565 fixU = true; 00566 if(fixU) { 00567 umin = lowu; 00568 umax = highu; 00569 min_du = umax-umin; 00570 best_startI = delta_startI; 00571 best_stopI = delta_stopI; 00572 } 00573 bool fixV = false; 00574 int delta_startJ = abs(startJ - (order_[1]+1)/2); 00575 int delta_stopJ = abs(stopJ - (order_[1]+1)/2); 00576 if( dv < min_dv ) 00577 fixV = true; 00578 if( dv == min_dv && delta_startJ <= best_startJ && delta_stopJ <= best_stopJ ) 00579 fixV = true; 00580 if(fixV) { 00581 vmin = lowv; 00582 vmax = highv; 00583 min_dv = vmax-vmin; 00584 best_startJ = delta_startJ; 00585 best_stopJ = delta_stopJ; 00586 } 00587 } 00588 00589 if(!only_insert_span_v_line) 00590 lines.push_back(new Meshline(true, (e->vmin() + e->vmax())/2.0, umin, umax, 1)); 00591 00592 if(!only_insert_span_u_line) 00593 lines.push_back(new Meshline(false, (e->umin() + e->umax())/2.0, vmin, vmax, 1)); 00594 00595 } 00596 00597 /************************************************************************************************************************/ 00605 void LRSplineSurface::getFullspanLines(int iEl, std::vector<Meshline*>& lines) { 00606 Element *e = element_[iEl]; 00607 double umin = e->umin(); 00608 double umax = e->umax(); 00609 double vmin = e->vmin(); 00610 double vmax = e->vmax(); 00611 bool only_insert_span_u_line = (vmax-vmin) >= maxAspectRatio_*(umax-umin); 00612 bool only_insert_span_v_line = (umax-umin) >= maxAspectRatio_*(vmax-vmin); 00613 // loop over all supported B-splines and make sure that everyone is covered by meshline 00614 for(Basisfunction *b : e->support()) { 00615 umin = (umin > (*b).getParmin(0)) ? (*b).getParmin(0) : umin; 00616 umax = (umax < (*b).getParmax(0)) ? (*b).getParmax(0) : umax; 00617 vmin = (vmin > (*b).getParmin(1)) ? (*b).getParmin(1) : vmin; 00618 vmax = (vmax < (*b).getParmax(1)) ? (*b).getParmax(1) : vmax; 00619 } 00620 if(!only_insert_span_v_line) 00621 lines.push_back(new Meshline(true, (e->vmin() + e->vmax())/2.0, umin, umax, 1)); 00622 00623 if(!only_insert_span_u_line) 00624 lines.push_back(new Meshline(false, (e->umin() + e->umax())/2.0, vmin, vmax, 1)); 00625 } 00626 00627 /************************************************************************************************************************/ 00635 void LRSplineSurface::getStructMeshLines(Basisfunction *b, std::vector<Meshline*>& lines) { 00636 double umin = b->getParmin(0); 00637 double umax = b->getParmax(0); 00638 double vmin = b->getParmin(1); 00639 double vmax = b->getParmax(1); 00640 00641 // find the largest knotspan in this function 00642 double max_du = 0; 00643 double max_dv = 0; 00644 for(int j=0; j<order_[0]; j++) { 00645 double du = (*b)[0][j+1]-(*b)[0][j]; 00646 bool isZeroSpan = fabs(du) < DOUBLE_TOL ; 00647 max_du = (isZeroSpan || max_du>du) ? max_du : du; 00648 } 00649 for(int j=0; j<order_[1]; j++) { 00650 double dv = (*b)[1][j+1]-(*b)[1][j]; 00651 bool isZeroSpan = fabs(dv) < DOUBLE_TOL ; 00652 max_dv = (isZeroSpan || max_dv>dv) ? max_dv : dv; 00653 } 00654 00655 // to keep as "square" basis function as possible, only insert 00656 // into the largest knot spans 00657 for(int j=0; j<order_[0]; j++) { 00658 double du = (*b)[0][j+1]-(*b)[0][j]; 00659 if( fabs(du-max_du) < DOUBLE_TOL ) 00660 lines.push_back(new Meshline(false, ((*b)[0][j] + (*b)[0][j+1])/2.0, vmin, vmax,1)); 00661 } 00662 for(int j=0; j<order_[1]; j++) { 00663 double dv = (*b)[1][j+1]-(*b)[1][j]; 00664 if( fabs(dv-max_dv) < DOUBLE_TOL ) 00665 lines.push_back(new Meshline(true, ((*b)[1][j] + (*b)[1][j+1])/2.0, umin, umax,1)); 00666 } 00667 } 00668 00669 /************************************************************************************************************************/ 00674 void LRSplineSurface::refineBasisFunction(int index) { 00675 std::vector<int> tmp = std::vector<int>(1, index); 00676 refineBasisFunction(tmp); 00677 } 00678 00679 /************************************************************************************************************************/ 00686 void LRSplineSurface::refineBasisFunction(const std::vector<int> &indices) { 00687 std::vector<int> sortedInd(indices); 00688 std::sort(sortedInd.begin(), sortedInd.end()); 00689 std::vector<Meshline*> newLines; 00690 00691 /* first retrieve all meshlines needed */ 00692 int ib = 0; 00693 HashSet_iterator<Basisfunction*> it = basis_.begin(); 00694 for(uint i=0; i<sortedInd.size(); i++) { 00695 while(ib < sortedInd[i]) { 00696 ++ib; 00697 ++it; 00698 } 00699 getStructMeshLines(*it,newLines); 00700 } 00701 00702 /* Do the actual refinement */ 00703 for(uint i=0; i<newLines.size(); i++) { 00704 Meshline *m = newLines[i]; 00705 insert_line(!m->is_spanning_u(), m->const_par_, m->start_, m->stop_, refKnotlineMult_); 00706 } 00707 00708 /* do a posteriori fixes to ensure a proper mesh */ 00709 aPosterioriFixes(); 00710 00711 /* exit cleanly be deleting all temporary new lines */ 00712 for(uint i=0; i<newLines.size(); i++) 00713 delete newLines[i]; 00714 } 00715 00716 /************************************************************************************************************************/ 00721 void LRSplineSurface::refineElement(int index) { 00722 std::vector<int> tmp = std::vector<int>(1, index); 00723 refineElement(tmp); 00724 } 00725 00726 /************************************************************************************************************************/ 00734 void LRSplineSurface::refineElement(const std::vector<int> &indices) { 00735 std::vector<Meshline*> newLines; 00736 00737 /* first retrieve all meshlines needed */ 00738 for(uint i=0; i<indices.size(); i++) { 00739 if(refStrat_ == LR_MINSPAN) 00740 getMinspanLines(indices[i],newLines); 00741 else 00742 getFullspanLines(indices[i],newLines); 00743 } 00744 00745 /* Do the actual refinement */ 00746 for(uint i=0; i<newLines.size(); i++) { 00747 Meshline *m = newLines[i]; 00748 insert_line(!m->is_spanning_u(), m->const_par_, m->start_, m->stop_, refKnotlineMult_); 00749 } 00750 00751 /* do a posteriori fixes to ensure a proper mesh */ 00752 aPosterioriFixes(); 00753 00754 /* exit cleanly be deleting all temporary new lines */ 00755 for(uint i=0; i<newLines.size(); i++) 00756 delete newLines[i]; 00757 } 00758 00759 void LRSplineSurface::refineByDimensionIncrease(const std::vector<double> &errPerElement, double beta) { 00760 Element *e; 00761 /* accumulate the error & index - vector */ 00762 std::vector<IndexDouble> errors; 00763 if(refStrat_ == LR_STRUCTURED_MESH) { // error per-function 00764 int i=0; 00765 for(Basisfunction *b : basis_) { 00766 errors.push_back(IndexDouble(0.0, i)); 00767 for(int j=0; j<b->nSupportedElements(); j++) { 00768 e = *(b->supportedElementBegin() + j); 00769 errors[i].first += errPerElement[e->getId()]; 00770 } 00771 i++; 00772 } 00773 } else { 00774 for(uint i=0; i<element_.size(); i++) 00775 errors.push_back(IndexDouble(errPerElement[i], i)); 00776 } 00777 00778 /* sort errors */ 00779 std::sort(errors.begin(), errors.end(), std::greater<IndexDouble>()); 00780 00781 /* first retrieve all possible meshlines needed */ 00782 std::vector<std::vector<Meshline*> > newLines(errors.size(), std::vector<Meshline*>(0)); 00783 for(uint i=0; i<errors.size(); i++) { 00784 if(refStrat_ == LR_MINSPAN) 00785 getMinspanLines(errors[i].second, newLines[i]); 00786 else if(refStrat_ == LR_FULLSPAN) 00787 getFullspanLines(errors[i].second, newLines[i]); 00788 else if(refStrat_ == LR_STRUCTURED_MESH) { 00789 Basisfunction *b = getBasisfunction(errors[i].second); 00790 getStructMeshLines(b, newLines[i]); 00791 } 00792 // note that this is an excessive loop as it computes the meshlines for ALL elements, 00793 // but we're only going to use a small part of this. 00794 } 00795 00796 /* Do the actual refinement */ 00797 int target_n_functions = ceil(basis_.size()*(1+beta)); 00798 int i=0; 00799 while( basis_.size() < target_n_functions ) { 00800 for(uint j=0; j<newLines[i].size(); j++) { 00801 Meshline *m = newLines[i][j]; 00802 insert_line(!m->is_spanning_u(), m->const_par_, m->start_, m->stop_, refKnotlineMult_); 00803 } 00804 i++; 00805 } 00806 00807 /* do a posteriori fixes to ensure a proper mesh */ 00808 aPosterioriFixes(); 00809 00810 /* exit cleanly be deleting all temporary new lines */ 00811 for(uint i=0; i<newLines.size(); i++) 00812 for(uint j=0; j<newLines[i].size(); j++) 00813 delete newLines[i][j]; 00814 } 00815 00816 void LRSplineSurface::aPosterioriFixes() { 00817 std::vector<Meshline*> *newLines = NULL; 00818 int nFunc; 00819 do { 00820 nFunc = basis_.size(); 00821 if(doCloseGaps_) 00822 this->closeGaps(newLines); 00823 if(maxTjoints_ > 0) 00824 this->enforceMaxTjoints(newLines); 00825 if(doAspectRatioFix_) 00826 this->enforceMaxAspectRatio(newLines); 00827 } while(nFunc != basis_.size()); 00828 } 00829 00830 00831 void LRSplineSurface::closeGaps(std::vector<Meshline*>* newLines) { 00832 std::vector<double> start_v; 00833 std::vector<double> stop_v ; 00834 std::vector<double> const_u ; 00835 00836 std::vector<double> start_u; 00837 std::vector<double> stop_u ; 00838 std::vector<double> const_v ; 00839 for(uint i=0; i<element_.size(); i++) { 00840 double umin = element_[i]->umin(); 00841 double umax = element_[i]->umax(); 00842 double vmin = element_[i]->vmin(); 00843 double vmax = element_[i]->vmax(); 00844 std::vector<double> left, right, top, bottom; 00845 for(uint j=0; j<meshline_.size(); j++) { 00846 Meshline *m = meshline_[j]; 00847 if( m->span_u_line_ && 00848 m->const_par_ > vmin && 00849 m->const_par_ < vmax) { 00850 if(m->start_ == umax) 00851 right.push_back(m->const_par_); 00852 else if(m->stop_ == umin) 00853 left.push_back(m->const_par_); 00854 } else if(!m->span_u_line_ && 00855 m->const_par_ > umin && 00856 m->const_par_ < umax) { 00857 if(m->start_ == vmax) 00858 top.push_back(m->const_par_); 00859 else if(m->stop_ == vmin) 00860 bottom.push_back(m->const_par_); 00861 } 00862 } 00863 for(uint j=0; j<left.size(); j++) 00864 for(uint k=0; k<right.size(); k++) 00865 if(left[j] == right[k]) { 00866 const_v.push_back(left[j]); 00867 start_u.push_back(umin); 00868 stop_u.push_back(umax); 00869 } 00870 for(uint j=0; j<top.size(); j++) 00871 for(uint k=0; k<bottom.size(); k++) 00872 if(top[j] == bottom[k]) { 00873 const_u.push_back(top[j]); 00874 start_v.push_back(vmin); 00875 stop_v.push_back(vmax); 00876 } 00877 } 00878 Meshline* m; 00879 for(uint i=0; i<const_u.size(); i++) { 00880 m = insert_const_u_edge(const_u[i], start_v[i], stop_v[i], refKnotlineMult_); 00881 if(newLines != NULL) 00882 newLines->push_back(m->copy()); 00883 } 00884 for(uint i=0; i<const_v.size(); i++) { 00885 m = insert_const_v_edge(const_v[i], start_u[i], stop_u[i], refKnotlineMult_); 00886 if(newLines != NULL) 00887 newLines->push_back(m->copy()); 00888 } 00889 } 00890 00891 void LRSplineSurface::enforceMaxTjoints(std::vector<Meshline*> *newLines) { 00892 bool someFix = true; 00893 while(someFix) { 00894 someFix = false; 00895 for(uint i=0; i<element_.size(); i++) { 00896 double umin = element_[i]->umin(); 00897 double umax = element_[i]->umax(); 00898 double vmin = element_[i]->vmin(); 00899 double vmax = element_[i]->vmax(); 00900 std::vector<double> left, right, top, bottom; 00901 for(uint j=0; j<meshline_.size(); j++) { 00902 Meshline *m = meshline_[j]; 00903 if( m->span_u_line_ && 00904 m->const_par_ > vmin && 00905 m->const_par_ < vmax) { 00906 if(m->start_ == umax) 00907 right.push_back(m->const_par_); 00908 else if(m->stop_ == umin) 00909 left.push_back(m->const_par_); 00910 } else if(!m->span_u_line_ && 00911 m->const_par_ > umin && 00912 m->const_par_ < umax) { 00913 if(m->start_ == vmax) 00914 top.push_back(m->const_par_); 00915 else if(m->stop_ == vmin) 00916 bottom.push_back(m->const_par_); 00917 } 00918 } 00919 Meshline *m; 00920 double best = DBL_MAX; 00921 int bi = -1; 00922 if(left.size() > (uint) maxTjoints_) { 00923 for(uint j=0; j<left.size(); j++) { 00924 if(fabs(left[j] - (vmin+vmax)/2) < best) { 00925 best = fabs(left[j] - (vmin+vmax)/2); 00926 bi = j; 00927 } 00928 } 00929 m = insert_const_v_edge(left[bi], umin, umax, refKnotlineMult_); 00930 if(newLines != NULL) 00931 newLines->push_back(m->copy()); 00932 someFix = true; 00933 continue; 00934 } 00935 if(right.size() > (uint) maxTjoints_) { 00936 for(uint j=0; j<right.size(); j++) { 00937 if(fabs(right[j] - (vmin+vmax)/2) < best) { 00938 best = fabs(right[j] - (vmin+vmax)/2); 00939 bi = j; 00940 } 00941 } 00942 m = insert_const_v_edge(right[bi], umin, umax, refKnotlineMult_); 00943 if(newLines != NULL) 00944 newLines->push_back(m->copy()); 00945 someFix = true; 00946 continue; 00947 } 00948 if(top.size() > (uint) maxTjoints_) { 00949 for(uint j=0; j<top.size(); j++) { 00950 if(fabs(top[j] - (umin+umax)/2) < best) { 00951 best = fabs(top[j] - (umin+umax)/2); 00952 bi = j; 00953 } 00954 } 00955 m = insert_const_u_edge(top[bi], vmin, vmax, refKnotlineMult_); 00956 if(newLines != NULL) 00957 newLines->push_back(m->copy()); 00958 someFix = true; 00959 continue; 00960 } 00961 if(bottom.size() > (uint) maxTjoints_) { 00962 for(uint j=0; j<bottom.size(); j++) { 00963 if(fabs(bottom[j] - (umin+umax)/2) < best) { 00964 best = fabs(bottom[j] - (umin+umax)/2); 00965 bi = j; 00966 } 00967 } 00968 m = insert_const_u_edge(bottom[bi], vmin, vmax, refKnotlineMult_); 00969 if(newLines != NULL) 00970 newLines->push_back(m->copy()); 00971 someFix = true; 00972 continue; 00973 } 00974 } 00975 } 00976 } 00977 00978 void LRSplineSurface::enforceMaxAspectRatio(std::vector<Meshline*>* newLines) { 00979 bool somethingFixed = true; 00980 while(somethingFixed) { 00981 somethingFixed = false; 00982 for(uint i=0; i<element_.size(); i++) { 00983 double umin = element_[i]->umin(); 00984 double umax = element_[i]->umax(); 00985 double vmin = element_[i]->vmin(); 00986 double vmax = element_[i]->vmax(); 00987 bool insert_const_u = umax-umin > maxAspectRatio_*(vmax-vmin); 00988 bool insert_const_v = vmax-vmin > maxAspectRatio_*(umax-umin); 00989 if( insert_const_u || insert_const_v ) { 00990 std::vector<Meshline*> splitLines; // should always contain exactly one meshline on function return 00991 if(refStrat_ == LR_MINSPAN) 00992 getMinspanLines(i, splitLines); 00993 else 00994 getFullspanLines(i, splitLines); 00995 00996 00997 Meshline *m, *msplit; 00998 msplit = splitLines.front(); 00999 01000 m = insert_line(!msplit->is_spanning_u(), msplit->const_par_, msplit->start_, msplit->stop_, refKnotlineMult_); 01001 if(newLines != NULL) 01002 newLines->push_back(m->copy()); 01003 01004 delete msplit; 01005 somethingFixed = true; 01006 break; 01007 } 01008 } 01009 } 01010 } 01011 01012 Meshline* LRSplineSurface::insert_const_u_edge(double u, double start_v, double stop_v, int multiplicity) { 01013 return insert_line(true, u, start_v, stop_v, multiplicity); 01014 } 01015 01016 Meshline* LRSplineSurface::insert_line(bool const_u, double const_par, double start, double stop, int multiplicity) { 01017 Meshline *newline; 01018 #ifdef TIME_LRSPLINE 01019 PROFILE("insert_line()"); 01020 #endif 01021 { // check if the line is an extension or a merging of existing lines 01022 #ifdef TIME_LRSPLINE 01023 PROFILE("line verification"); 01024 #endif 01025 newline = new Meshline(!const_u, const_par, start, stop, multiplicity); 01026 newline->type_ = NEWLINE; 01027 for(uint i=0; i<meshline_.size(); i++) { 01028 // if newline overlaps any existing ones (may be multiple existing ones) 01029 // let newline be the entire length of all merged and delete the unused ones 01030 01031 if(meshline_[i]->is_spanning_u() != const_u && fabs(meshline_[i]->const_par_-const_par)<DOUBLE_TOL && 01032 meshline_[i]->start_ <= stop && meshline_[i]->stop_ >= start) { // meshline_[i] overlaps with newline 01033 01034 if(meshline_[i]->start_ <= start && 01035 meshline_[i]->stop_ >= stop ) { // newline completely contained in meshline_[i] 01036 01037 if(meshline_[i]->multiplicity_ < newline->multiplicity_) { // increasing multiplicity 01038 if(meshline_[i]->start_ == start && 01039 meshline_[i]->stop_ == stop ) { // increasing the mult of the entire line 01040 01041 // keeping newline, getting rid of the old line 01042 delete meshline_[i]; 01043 meshline_.erase(meshline_.begin() + i); 01044 i--; 01045 01046 } else { // increasing multiplicity of partial line 01047 // do nothing. Keep the entire length meshline_[i], and add newline 01048 01049 } 01050 01051 } else { // line exist already, do nothing 01052 delete newline; 01053 return meshline_[i]; 01054 } 01055 } else { // newline overlaps meshline_[i]. Keep (and update) newline, delete meshline_[i] 01056 01057 // update refinement type (for later analysis of linear independence) 01058 if(newline->type_ == ELONGATION) // overlaps two existing lines => MERGING 01059 newline->type_ = MERGING; 01060 else if(newline->type_ != MERGING) // overlaps one existing line => ELONGATION 01061 newline->type_ = ELONGATION; 01062 01063 // update the length of the line with the lowest multiplicity 01064 if(meshline_[i]->multiplicity_ < newline->multiplicity_) { 01065 if(meshline_[i]->start_ > start) meshline_[i]->start_ = newline->start_; 01066 if(meshline_[i]->stop_ < stop ) meshline_[i]->stop_ = newline->stop_; 01067 01068 } else if(meshline_[i]->multiplicity_ > newline->multiplicity_) { 01069 if(meshline_[i]->start_ < start) newline->start_ = meshline_[i]->start_; 01070 if(meshline_[i]->stop_ > stop ) newline->stop_ = meshline_[i]->stop_; 01071 01072 } else { // for equal mult, we only keep newline and remove the previous line 01073 if(meshline_[i]->start_ < start) newline->start_ = meshline_[i]->start_; 01074 if(meshline_[i]->stop_ > stop ) newline->stop_ = meshline_[i]->stop_; 01075 01076 // keeping newline, getting rid of the old line 01077 delete meshline_[i]; 01078 meshline_.erase(meshline_.begin() + i); 01079 i--; 01080 } 01081 01082 } 01083 } 01084 } 01085 } 01086 01087 HashSet<Basisfunction*> newFuncStp1, newFuncStp2; 01088 HashSet<Basisfunction*> removeFunc; 01089 01090 { // STEP 1: test EVERY function against the NEW meshline 01091 #ifdef TIME_LRSPLINE 01092 PROFILE("STEP 1"); 01093 #endif 01094 { 01095 #ifdef TIME_LRSPLINE 01096 PROFILE("S1-basissplit"); 01097 #endif 01098 for(Basisfunction* b : basis_) { 01099 if(newline->splits(b)) { 01100 int nKnots; 01101 if( ((nKnots=newline->nKnotsIn(b)) != newline->multiplicity_) ) { 01102 removeFunc.insert(b); 01103 split( const_u, b, const_par, newline->multiplicity_-nKnots, newFuncStp1 ); 01104 } 01105 } 01106 } 01107 for(Basisfunction* b : removeFunc) { 01108 basis_.erase(b); 01109 delete b; 01110 } 01111 } // end profiler 01112 { 01113 #ifdef TIME_LRSPLINE 01114 PROFILE("S1-elementsplit"); 01115 #endif 01116 for(uint i=0; i<element_.size(); i++) { 01117 if(newline->splits(element_[i])) 01118 element_.push_back(element_[i]->split(newline->is_spanning_u(), newline->const_par_)); 01119 } 01120 } // end profiler (elementsplit) 01121 } // end profiler (step 1) 01122 01123 { // STEP 2: test every NEW function against ALL old meshlines 01124 #ifdef TIME_LRSPLINE 01125 PROFILE("STEP 2"); 01126 #endif 01127 meshline_.push_back(newline); 01128 while(newFuncStp1.size() > 0) { 01129 Basisfunction *b = newFuncStp1.pop(); 01130 bool splitMore = false; 01131 for(Meshline *m : meshline_) { 01132 if(m->splits(b)) { 01133 int nKnots = m->nKnotsIn(b); 01134 if( nKnots != m->multiplicity_ ) { 01135 splitMore = true; 01136 split( !m->is_spanning_u(), b, m->const_par_, m->multiplicity_-nKnots, newFuncStp1); 01137 delete b; 01138 break; 01139 } 01140 } 01141 } 01142 if(!splitMore) 01143 basis_.insert(b); 01144 } 01145 } // end profiler (step 2) 01146 01147 return newline; 01148 } 01149 01150 Meshline* LRSplineSurface::insert_const_v_edge(double v, double start_u, double stop_u, int multiplicity) { 01151 return insert_line(false, v, start_u, stop_u, multiplicity); 01152 } 01153 01154 void LRSplineSurface::split(bool insert_in_u, Basisfunction* b, double new_knot, int multiplicity, HashSet<Basisfunction*> &newFunctions) { 01155 #ifdef TIME_LRSPLINE 01156 PROFILE("split()"); 01157 #endif 01158 01159 // create the new functions b1 and b2 01160 Basisfunction *b1, *b2; 01161 std::vector<double> knot = (insert_in_u) ? (*b)[0] : (*b)[1]; 01162 int p = (insert_in_u) ? order_[0] : order_[1]; 01163 int insert_index = 0; 01164 if(new_knot < knot[0] || knot[p] < new_knot) 01165 return ; 01166 while(knot[insert_index] < new_knot) 01167 insert_index++; 01168 double alpha1 = (insert_index == p) ? 1.0 : (new_knot-knot[0])/(knot[p-1]-knot[0]); 01169 double alpha2 = (insert_index == 1 ) ? 1.0 : (knot[p]-new_knot)/(knot[p]-knot[1]); 01170 double newKnot[p+2]; 01171 std::copy(knot.begin(), knot.end(), newKnot+1); 01172 newKnot[0] = new_knot; 01173 std::sort(newKnot, newKnot + p+2); 01174 if(insert_in_u) { 01175 b1 = new Basisfunction(newKnot , (*b)[1].begin(), b->cp(), b->dim(), order_[0], order_[1], b->w()*alpha1); 01176 b2 = new Basisfunction(newKnot+1, (*b)[1].begin(), b->cp(), b->dim(), order_[0], order_[1], b->w()*alpha2); 01177 } else { // insert in v 01178 b1 = new Basisfunction((*b)[0].begin(), newKnot , b->cp(), b->dim(), order_[0], order_[1], b->w()*alpha1); 01179 b2 = new Basisfunction((*b)[0].begin(), newKnot+1, b->cp(), b->dim(), order_[0], order_[1], b->w()*alpha2); 01180 } 01181 01182 // add any brand new functions and detect their support elements 01183 HashSet_iterator<Basisfunction*> it = basis_.find(b1); 01184 if(it != basis_.end()) { 01185 **it += *b1; 01186 delete b1; 01187 } else { 01188 it = newFunctions.find(b1); 01189 if(it != newFunctions.end()) { 01190 **it += *b1; 01191 delete b1; 01192 } else { 01193 updateSupport(b1, b->supportedElementBegin(), b->supportedElementEnd()); 01194 bool recursive_split = (multiplicity > 1) && ( ( insert_in_u && (*b1)[0][order_[0]]!=new_knot) || 01195 (!insert_in_u && (*b1)[1][order_[1]]!=new_knot) ); 01196 if(recursive_split) { 01197 split( insert_in_u, b1, new_knot, multiplicity-1, newFunctions); 01198 delete b1; 01199 } else { 01200 newFunctions.insert(b1); 01201 } 01202 } 01203 } 01204 it = basis_.find(b2); 01205 if(it != basis_.end()) { 01206 **it += *b2; 01207 delete b2; 01208 } else { 01209 it = newFunctions.find(b2); 01210 if(it != newFunctions.end()) { 01211 **it += *b2; 01212 delete b2; 01213 } else { 01214 updateSupport(b2, b->supportedElementBegin(), b->supportedElementEnd()); 01215 bool recursive_split = (multiplicity > 1) && ( ( insert_in_u && (*b2)[0][0]!=new_knot) || 01216 (!insert_in_u && (*b2)[1][0]!=new_knot) ); 01217 if(recursive_split) { 01218 split( insert_in_u, b2, new_knot, multiplicity-1, newFunctions); 01219 delete b2; 01220 } else { 01221 newFunctions.insert(b2); 01222 } 01223 } 01224 } 01225 01226 } 01227 01228 void LRSplineSurface::getGlobalKnotVector(std::vector<double> &knot_u, std::vector<double> &knot_v) const { 01229 getGlobalUniqueKnotVector(knot_u, knot_v); 01230 01231 // add in duplicates where apropriate 01232 for(uint i=0; i<knot_u.size(); i++) { 01233 for(uint j=0; j<meshline_.size(); j++) { 01234 if(!meshline_[j]->is_spanning_u() && meshline_[j]->const_par_==knot_u[i]) { 01235 for(int m=1; m<meshline_[j]->multiplicity_; m++) { 01236 knot_u.insert(knot_u.begin()+i, knot_u[i]); 01237 i++; 01238 } 01239 break; 01240 } 01241 } 01242 } 01243 for(uint i=0; i<knot_v.size(); i++) { 01244 for(uint j=0; j<meshline_.size(); j++) { 01245 if(meshline_[j]->is_spanning_u() && meshline_[j]->const_par_==knot_v[i]) { 01246 for(int m=1; m<meshline_[j]->multiplicity_; m++) { 01247 knot_v.insert(knot_v.begin()+i, knot_v[i]); 01248 i++; 01249 } 01250 break; 01251 } 01252 } 01253 } 01254 } 01255 01256 void LRSplineSurface::getGlobalUniqueKnotVector(std::vector<double> &knot_u, std::vector<double> &knot_v) const { 01257 knot_u.clear(); 01258 knot_v.clear(); 01259 // create a huge list of all line instances 01260 for(uint i=0; i<meshline_.size(); i++) { 01261 if(meshline_[i]->is_spanning_u()) 01262 knot_v.push_back(meshline_[i]->const_par_); 01263 else 01264 knot_u.push_back(meshline_[i]->const_par_); 01265 } 01266 01267 // sort them and remove duplicates 01268 std::vector<double>::iterator it; 01269 std::sort(knot_u.begin(), knot_u.end()); 01270 std::sort(knot_v.begin(), knot_v.end()); 01271 it = std::unique(knot_u.begin(), knot_u.end()); 01272 knot_u.resize(it-knot_u.begin()); 01273 it = std::unique(knot_v.begin(), knot_v.end()); 01274 knot_v.resize(it-knot_v.begin()); 01275 } 01276 01277 void LRSplineSurface::updateSupport(Basisfunction *f, 01278 std::vector<Element*>::iterator start, 01279 std::vector<Element*>::iterator end ) { 01280 #ifdef TIME_LRSPLINE 01281 PROFILE("update support"); 01282 #endif 01283 std::vector<Element*>::iterator it; 01284 for(it=start; it!=end; it++) 01285 if(f->addSupport(*it)) // this tests for overlapping as well as updating 01286 (*it)->addSupportFunction(f); 01287 } 01288 01289 void LRSplineSurface::updateSupport(Basisfunction *f) { 01290 updateSupport(f, element_.begin(), element_.end()); 01291 } 01292 01293 bool LRSplineSurface::isLinearIndepByOverloading(bool verbose) { 01294 std::vector<Basisfunction*> overloaded; 01295 std::vector<int> singleElms; 01296 std::vector<int> multipleElms; 01297 std::vector<int> singleBasis; 01298 std::vector<int> multipleBasis; 01299 // reset overload counts and initialize overloaded basisfunction set 01300 for(uint i=0; i<element_.size(); i++) 01301 element_[i]->resetOverloadCount(); 01302 for(Basisfunction *b : basis_) 01303 if(b->isOverloaded()) 01304 overloaded.push_back(b); 01305 01306 int lastOverloadCount = overloaded.size(); 01307 int iterationCount = 0 ; 01308 do { 01309 lastOverloadCount = overloaded.size(); 01310 // reset overload counts and plotting vectors 01311 for(uint i=0; i<element_.size(); i++) 01312 element_[i]->resetOverloadCount(); 01313 singleBasis.clear(); 01314 multipleBasis.clear(); 01315 singleElms.clear(); 01316 multipleElms.clear(); 01317 for(uint i=0; i<overloaded.size(); i++) 01318 for(Element* e : overloaded[i]->support()) 01319 e->incrementOverloadCount(); 01320 for(uint i=0; i<element_.size(); i++) { 01321 if(element_[i]->getOverloadCount() > 1) 01322 multipleElms.push_back(i); 01323 if(element_[i]->getOverloadCount() == 1) 01324 singleElms.push_back(i); 01325 } 01326 for(uint i=0; i<overloaded.size(); i++) { 01327 if(overloaded[i]->getOverloadCount() > 1) 01328 multipleBasis.push_back(overloaded[i]->getId()); 01329 if(overloaded[i]->getOverloadCount() == 1) { 01330 singleBasis.push_back(overloaded[i]->getId()); 01331 overloaded.erase(overloaded.begin() + i); 01332 i--; 01333 } 01334 } 01335 int nOverloadedElms = singleElms.size() + multipleElms.size(); 01336 if(verbose) { 01337 std::cout << "Iteration #"<< iterationCount << "\n"; 01338 std::cout << "Overloaded elements : " << nOverloadedElms 01339 << " ( " << singleElms.size() << " + " 01340 << multipleElms.size() << ")\n"; 01341 std::cout << "Overloaded B-splines : " << lastOverloadCount 01342 << " ( " << overloaded.size() << " + " 01343 << (lastOverloadCount-overloaded.size()) << ")\n"; 01344 std::cout << "-----------------------------------------------\n\n"; 01345 01346 char filename[256]; 01347 std::ofstream out; 01348 sprintf(filename, "elements_%03d.eps", iterationCount); 01349 out.open(filename); 01350 setElementColor(0.6, 0.6, 0.6); 01351 writePostscriptElements(out, 2,2,false, &singleElms); 01352 setElementColor(0.9, 0.3, 0.15); 01353 writePostscriptElements(out, 2,2,true, &multipleElms); 01354 out.close(); 01355 01356 sprintf(filename, "functions_%03d.eps", iterationCount); 01357 out.open(filename); 01358 setBasisColor(0.6, 0.6, 0.6); 01359 setSelectedBasisColor(0.9, 0.83, 0.05); 01360 writePostscriptFunctionSpace(out, &singleBasis, true, false); 01361 setSelectedBasisColor(1.0, 0.10, 0.00); 01362 writePostscriptFunctionSpace(out, &multipleBasis, false, true); 01363 out.close(); 01364 } 01365 01366 01367 iterationCount++; 01368 } while((uint) lastOverloadCount != overloaded.size()); 01369 01370 return overloaded.size() == 0; 01371 } 01372 01373 bool LRSplineSurface::isLinearIndepByMappingMatrix(bool verbose) const { 01374 #ifdef TIME_LRSPLINE 01375 PROFILE("Linear independent)"); 01376 #endif 01377 // try and figure out this thing by the projection matrix C 01378 01379 std::vector<double> knots_u, knots_v; 01380 getGlobalKnotVector(knots_u, knots_v); 01381 int nmb_bas = basis_.size(); 01382 int n1 = knots_u.size() - order_[0]; 01383 int n2 = knots_v.size() - order_[1]; 01384 int fullDim = n1*n2; 01385 bool fullVerbose = fullDim < 30 && nmb_bas < 50; 01386 bool sparseVerbose = fullDim < 250 && nmb_bas < 100; 01387 01388 std::vector<std::vector<boost::rational<long long> > > C; // rational projection matrix 01389 01390 // scaling factor to ensure that all knots are integers (assuming all multiplum of smallest knot span) 01391 double smallKnotU = DBL_MAX; 01392 double smallKnotV = DBL_MAX; 01393 for(uint i=0; i<knots_u.size()-1; i++) 01394 if(knots_u[i+1]-knots_u[i] < smallKnotU && knots_u[i+1] != knots_u[i]) 01395 smallKnotU = knots_u[i+1]-knots_u[i]; 01396 for(uint i=0; i<knots_v.size()-1; i++) 01397 if(knots_v[i+1]-knots_v[i] < smallKnotV && knots_v[i+1] != knots_v[i]) 01398 smallKnotV = knots_v[i+1]-knots_v[i]; 01399 01400 // HashSet_iterator<Basisfunction*> it_begin = basis_.begin(); 01401 // HashSet_iterator<Basisfunction*> it_end = basis_.end(); 01402 HashSet_const_iterator<Basisfunction*> cit_begin = basis_.begin(); 01403 HashSet_const_iterator<Basisfunction*> cit_end = basis_.end(); 01404 HashSet_const_iterator<Basisfunction*> it; 01405 // for(Basisfunction *b : basis_) { 01406 for(it=cit_begin; it != cit_end; ++it) { 01407 Basisfunction *b = *it; 01408 int startU, startV; 01409 std::vector<double> locKnotU((*b)[0].begin(), (*b)[0].end()); 01410 std::vector<double> locKnotV((*b)[1].begin(), (*b)[1].end()); 01411 01412 for(startU=knots_u.size(); startU-->0; ) 01413 if(knots_u[startU] == (*b)[0][0]) break; 01414 for(int j=0; j<order_[0]; j++) { 01415 if(knots_u[startU] == (*b)[0][j]) startU--; 01416 else break; 01417 } 01418 startU++; 01419 for(startV=knots_v.size(); startV-->0; ) 01420 if(knots_v[startV] == (*b)[1][0]) break; 01421 for(int j=0; j<order_[1]; j++) { 01422 if(knots_v[startV] == (*b)[1][j]) startV--; 01423 else break; 01424 } 01425 startV++; 01426 01427 std::vector<boost::rational<long long> > rowU(1,1), rowV(1,1); 01428 int curU = startU+1; 01429 for(uint j=0; j<locKnotU.size()-1; j++, curU++) { 01430 if(locKnotU[j+1] != knots_u[curU]) { 01431 std::vector<boost::rational<long long> > newRowU(rowU.size()+1, boost::rational<long long>(0)); 01432 for(uint k=0; k<rowU.size(); k++) { 01433 #define U(x) ((long long) (locKnotU[x+k]/smallKnotU + 0.5)) 01434 long long z = (long long) (knots_u[curU] / smallKnotU + 0.5); 01435 int p = order_[0]-1; 01436 if(z < U(0) || z > U(p+1)) { 01437 newRowU[k] = rowU[k]; 01438 continue; 01439 } 01440 boost::rational<long long> alpha1 = (U(p) <= z ) ? 1 : boost::rational<long long>( z - U(0), U(p) - U(0)); 01441 boost::rational<long long> alpha2 = (z <= U(1)) ? 1 : boost::rational<long long>( U(p+1) - z , U(p+1) - U(1)); 01442 newRowU[k] += rowU[k]*alpha1; 01443 newRowU[k+1] += rowU[k]*alpha2; 01444 #undef U 01445 } 01446 locKnotU.insert(locKnotU.begin()+(j+1), knots_u[curU]); 01447 rowU = newRowU; 01448 } 01449 } 01450 int curV = startV+1; 01451 for(uint j=0; j<locKnotV.size()-1; j++, curV++) { 01452 if(locKnotV[j+1] != knots_v[curV]) { 01453 std::vector<boost::rational<long long> > newRowV(rowV.size()+1, boost::rational<long long>(0)); 01454 for(uint k=0; k<rowV.size(); k++) { 01455 #define V(x) ((long long) (locKnotV[x+k]/smallKnotV + 0.5)) 01456 long long z = (long long) (knots_v[curV] / smallKnotV + 0.5); 01457 int p = order_[1]-1; 01458 if(z < V(0) || z > V(p+1)) { 01459 newRowV[k] = rowV[k]; 01460 continue; 01461 } 01462 boost::rational<long long> alpha1 = (V(p) <= z ) ? 1 : boost::rational<long long>( z - V(0), V(p) - V(0)); 01463 boost::rational<long long> alpha2 = (z <= V(1)) ? 1 : boost::rational<long long>( V(p+1) - z , V(p+1) - V(1)); 01464 newRowV[k] += rowV[k]*alpha1; 01465 newRowV[k+1] += rowV[k]*alpha2; 01466 #undef V 01467 } 01468 locKnotV.insert(locKnotV.begin()+(j+1), knots_v[curV]); 01469 rowV= newRowV; 01470 } 01471 } 01472 std::vector<boost::rational<long long> > totalRow(fullDim, boost::rational<long long>(0)); 01473 for(uint i1=0; i1<rowU.size(); i1++) 01474 for(uint i2=0; i2<rowV.size(); i2++) 01475 totalRow[(startV+i2)*n1 + (startU+i1)] = rowV[i2]*rowU[i1]; 01476 01477 C.push_back(totalRow); 01478 } 01479 01480 if(verbose && sparseVerbose) { 01481 for(uint i=0; i<C.size(); i++) { 01482 std::cout << "|"; 01483 for(uint j=0; j<C[i].size(); j++) 01484 if(C[i][j]==0) { 01485 if(fullVerbose) 01486 std::cout << "\t"; 01487 else if(sparseVerbose) 01488 std::cout << " "; 01489 } else { 01490 if(fullVerbose) 01491 std::cout << C[i][j] << "\t"; 01492 else if(sparseVerbose) 01493 std::cout << "x"; 01494 } 01495 std::cout << "|\n"; 01496 } 01497 std::cout << std::endl; 01498 } 01499 01500 // gauss-jordan elimination 01501 int zeroColumns = 0; 01502 for(uint i=0; i<C.size() && i+zeroColumns<C[i].size(); i++) { 01503 boost::rational<long long> maxPivot = 0; 01504 int maxI = -1; 01505 for(uint j=i; j<C.size(); j++) { 01506 if(abs(C[j][i+zeroColumns]) > maxPivot) { 01507 maxPivot = abs(C[j][i+zeroColumns]); 01508 maxI = j; 01509 } 01510 } 01511 if(maxI == -1) { 01512 i--; 01513 zeroColumns++; 01514 continue; 01515 } 01516 std::vector<boost::rational<long long> > tmp = C[i]; 01517 C[i] = C[maxI]; 01518 C[maxI] = tmp; 01519 for(uint j=i+1; j<C.size(); j++) { 01520 // if(j==i) continue; 01521 boost::rational<long long> scale = C[j][i+zeroColumns] / C[i][i+zeroColumns]; 01522 if(scale != 0) { 01523 // std::cout << scale << " x row(" << i << ") added to row " << j << std::endl; 01524 for(uint k=i+zeroColumns; k<C[j].size(); k++) 01525 C[j][k] -= C[i][k] * scale; 01526 } 01527 } 01528 } 01529 01530 if(verbose && sparseVerbose) { 01531 for(uint i=0; i<C.size(); i++) { 01532 std::cout << "|"; 01533 for(uint j=0; j<C[i].size(); j++) 01534 if(C[i][j]==0) { 01535 if(fullVerbose) 01536 std::cout << "\t"; 01537 else if(sparseVerbose) 01538 std::cout << " "; 01539 } else { 01540 if(fullVerbose) 01541 std::cout << C[i][j] << "\t"; 01542 else if(sparseVerbose) 01543 std::cout << "x"; 01544 } 01545 std::cout << "|\n"; 01546 } 01547 std::cout << std::endl; 01548 } 01549 01550 int rank = (nmb_bas < n1*n2-zeroColumns) ? nmb_bas : n1*n2-zeroColumns; 01551 if(verbose) { 01552 std::cout << "Matrix size : " << nmb_bas << " x " << n1*n2 << std::endl; 01553 std::cout << "Matrix rank : " << rank << std::endl; 01554 } 01555 01556 return rank == nmb_bas; 01557 } 01558 01559 void LRSplineSurface::getNullSpace(std::vector<std::vector<boost::rational<long long> > >& nullspace) const { 01560 #ifdef TIME_LRSPLINE 01561 PROFILE("Linear independent)"); 01562 #endif 01563 // try and figure out this thing by the projection matrix C 01564 01565 std::vector<double> knots_u, knots_v; 01566 getGlobalKnotVector(knots_u, knots_v); 01567 int nmb_bas = basis_.size(); 01568 int n1 = knots_u.size() - order_[0]; 01569 int n2 = knots_v.size() - order_[1]; 01570 int fullDim = n1*n2; 01571 01572 std::vector<std::vector<boost::rational<long long> > > Ct; // rational projection matrix (transpose of this) 01573 01574 // scaling factor to ensure that all knots are integers (assuming all multiplum of smallest knot span) 01575 double smallKnotU = DBL_MAX; 01576 double smallKnotV = DBL_MAX; 01577 for(uint i=0; i<knots_u.size()-1; i++) 01578 if(knots_u[i+1]-knots_u[i] < smallKnotU && knots_u[i+1] != knots_u[i]) 01579 smallKnotU = knots_u[i+1]-knots_u[i]; 01580 for(uint i=0; i<knots_v.size()-1; i++) 01581 if(knots_v[i+1]-knots_v[i] < smallKnotV && knots_v[i+1] != knots_v[i]) 01582 smallKnotV = knots_v[i+1]-knots_v[i]; 01583 01584 // initialize C transpose with zeros 01585 for(int i=0; i<fullDim; i++) 01586 Ct.push_back(std::vector<boost::rational<long long> >(nmb_bas, boost::rational<long long>(0))); 01587 01588 int basCount = 0; 01589 for(Basisfunction *b : basis_) { 01590 int startU, startV; 01591 std::vector<double> locKnotU((*b)[0].begin(), (*b)[0].end()); 01592 std::vector<double> locKnotV((*b)[1].begin(), (*b)[1].end()); 01593 01594 for(startU=knots_u.size(); startU-->0; ) 01595 if(knots_u[startU] == (*(b))[0][0]) break; 01596 for(int j=0; j<order_[0]; j++) { 01597 if(knots_u[startU] == (*(b))[0][j]) startU--; 01598 else break; 01599 } 01600 startU++; 01601 for(startV=knots_v.size(); startV-->0; ) 01602 if(knots_v[startV] == (*(b))[1][0]) break; 01603 for(int j=0; j<order_[1]; j++) { 01604 if(knots_v[startV] == (*(b))[1][j]) startV--; 01605 else break; 01606 } 01607 startV++; 01608 01609 std::vector<boost::rational<long long> > rowU(1,1), rowV(1,1); 01610 int curU = startU+1; 01611 for(uint j=0; j<locKnotU.size()-1; j++, curU++) { 01612 if(locKnotU[j+1] != knots_u[curU]) { 01613 std::vector<boost::rational<long long> > newRowU(rowU.size()+1, boost::rational<long long>(0)); 01614 for(uint k=0; k<rowU.size(); k++) { 01615 #define U(x) ((long long) (locKnotU[x+k]/smallKnotU + 0.5)) 01616 long long z = (long long) (knots_u[curU] / smallKnotU + 0.5); 01617 int p = order_[0]-1; 01618 if(z < U(0) || z > U(p+1)) { 01619 newRowU[k] = rowU[k]; 01620 continue; 01621 } 01622 boost::rational<long long> alpha1 = (U(p) <= z ) ? 1 : boost::rational<long long>( z - U(0), U(p) - U(0)); 01623 boost::rational<long long> alpha2 = (z <= U(1)) ? 1 : boost::rational<long long>( U(p+1) - z , U(p+1) - U(1)); 01624 newRowU[k] += rowU[k]*alpha1; 01625 newRowU[k+1] += rowU[k]*alpha2; 01626 #undef U 01627 } 01628 locKnotU.insert(locKnotU.begin()+(j+1), knots_u[curU]); 01629 rowU = newRowU; 01630 } 01631 } 01632 int curV = startV+1; 01633 for(uint j=0; j<locKnotV.size()-1; j++, curV++) { 01634 if(locKnotV[j+1] != knots_v[curV]) { 01635 std::vector<boost::rational<long long> > newRowV(rowV.size()+1, boost::rational<long long>(0)); 01636 for(uint k=0; k<rowV.size(); k++) { 01637 #define V(x) ((long long) (locKnotV[x+k]/smallKnotV + 0.5)) 01638 long long z = (long long) (knots_v[curV] / smallKnotV + 0.5); 01639 int p = order_[1]-1; 01640 if(z < V(0) || z > V(p+1)) { 01641 newRowV[k] = rowV[k]; 01642 continue; 01643 } 01644 boost::rational<long long> alpha1 = (V(p) <= z ) ? 1 : boost::rational<long long>( z - V(0), V(p) - V(0)); 01645 boost::rational<long long> alpha2 = (z <= V(1)) ? 1 : boost::rational<long long>( V(p+1) - z , V(p+1) - V(1)); 01646 newRowV[k] += rowV[k]*alpha1; 01647 newRowV[k+1] += rowV[k]*alpha2; 01648 #undef V 01649 } 01650 locKnotV.insert(locKnotV.begin()+(j+1), knots_v[curV]); 01651 rowV= newRowV; 01652 } 01653 } 01654 for(uint i1=0; i1<rowU.size(); i1++) 01655 for(uint i2=0; i2<rowV.size(); i2++) 01656 Ct[(startV+i2)*n1 + (startU+i1)][basCount] = rowV[i2]*rowU[i1]; 01657 basCount++; 01658 } 01659 01660 01661 std::ofstream out; 01662 out.open("Ct.m"); 01663 out << "A = ["; 01664 for(uint i=0; i<Ct.size(); i++) { 01665 for(uint j=0; j<Ct[i].size(); j++) 01666 if(Ct[i][j] > 0) 01667 out << i << " " << j << " " << Ct[i][j] << ";\n"; 01668 } 01669 out << "];" << std::endl; 01670 out.close(); 01671 // gauss-jordan elimination 01672 int zeroColumns = 0; 01673 for(uint i=0; i<Ct.size() && i+zeroColumns<Ct[i].size(); i++) { 01674 boost::rational<long long> maxPivot(0); 01675 int maxI = -1; 01676 for(uint j=i; j<Ct.size(); j++) { 01677 if(abs(Ct[j][i+zeroColumns]) > maxPivot) { 01678 maxPivot = abs(Ct[j][i+zeroColumns]); 01679 maxI = j; 01680 } 01681 } 01682 if(maxI == -1) { 01683 i--; 01684 zeroColumns++; 01685 continue; 01686 } 01687 std::vector<boost::rational<long long> > tmp = Ct[i]; 01688 Ct[i] = Ct[maxI]; 01689 Ct[maxI] = tmp; 01690 boost::rational<long long> leading = Ct[i][i+zeroColumns]; 01691 for(uint j=i+zeroColumns; j<Ct[i].size(); j++) 01692 Ct[i][j] /= leading; 01693 01694 for(uint j=0; j<Ct.size(); j++) { 01695 if(i==j) continue; 01696 boost::rational<long long> scale = Ct[j][i+zeroColumns] / Ct[i][i+zeroColumns]; 01697 if(scale != 0) { 01698 for(uint k=i+zeroColumns; k<Ct[j].size(); k++) 01699 Ct[j][k] -= Ct[i][k] * scale; 01700 } 01701 } 01702 } 01703 01704 // figure out the null space 01705 nullspace.clear(); 01706 for(int i=0; i<zeroColumns; i++) { 01707 std::vector<boost::rational<long long> > newRow(nmb_bas, boost::rational<long long>(0)); 01708 newRow[nmb_bas-(i+1)] = 1; 01709 nullspace.push_back(newRow); 01710 } 01711 for(int i=0; i<zeroColumns; i++) 01712 for(int j=0; j<nmb_bas-zeroColumns; j++) 01713 nullspace[i][j] = -Ct[j][nmb_bas-(i+1)]; // no need to divide by C[j][j] since it's 1 01714 01715 } 01716 01717 bool LRSplineSurface::isLinearIndepByFloatingPointMappingMatrix(bool verbose) const { 01718 #ifdef TIME_LRSPLINE 01719 PROFILE("Linear independent)"); 01720 #endif 01721 // try and figure out this thing by the projection matrix C 01722 01723 std::vector<double> knots_u, knots_v; 01724 getGlobalKnotVector(knots_u, knots_v); 01725 int nmb_bas = basis_.size(); 01726 int n1 = knots_u.size() - order_[0]; 01727 int n2 = knots_v.size() - order_[1]; 01728 int fullDim = n1*n2; 01729 bool fullVerbose = fullDim < 30 && nmb_bas < 50; 01730 bool sparseVerbose = fullDim < 250 && nmb_bas < 100; 01731 01732 std::vector<std::vector<double > > C; // rational projection matrix 01733 01734 for(Basisfunction *b : basis_) { 01735 int startU, startV; 01736 std::vector<double> locKnotU((*b)[0].begin(), (*b)[0].end()); 01737 std::vector<double> locKnotV((*b)[1].begin(), (*b)[1].end()); 01738 01739 for(startU=knots_u.size(); startU-->0; ) 01740 if(knots_u[startU] == (*b)[0][0]) break; 01741 for(int j=0; j<order_[0]; j++) { 01742 if(knots_u[startU] == (*b)[0][j]) startU--; 01743 else break; 01744 } 01745 startU++; 01746 for(startV=knots_v.size(); startV-->0; ) 01747 if(knots_v[startV] == (*b)[1][0]) break; 01748 for(int j=0; j<order_[1]; j++) { 01749 if(knots_v[startV] == (*b)[1][j]) startV--; 01750 else break; 01751 } 01752 startV++; 01753 01754 std::vector<double > rowU(1,1), rowV(1,1); 01755 int curU = startU+1; 01756 for(uint j=0; j<locKnotU.size()-1; j++, curU++) { 01757 if(locKnotU[j+1] != knots_u[curU]) { 01758 std::vector<double > newRowU(rowU.size()+1, 0); 01759 for(uint k=0; k<rowU.size(); k++) { 01760 #define U(x) (locKnotU[x+k]) 01761 double z = knots_u[curU] ; 01762 int p = order_[0]-1; 01763 if(z < U(0) || z > U(p+1)) { 01764 newRowU[k] = rowU[k]; 01765 continue; 01766 } 01767 double alpha1 = (U(p) <= z ) ? 1 : ( z - U(0))/( U(p) - U(0)); 01768 double alpha2 = (z <= U(1)) ? 1 : ( U(p+1) - z )/( U(p+1) - U(1)); 01769 newRowU[k] += rowU[k]*alpha1; 01770 newRowU[k+1] += rowU[k]*alpha2; 01771 #undef U 01772 } 01773 locKnotU.insert(locKnotU.begin()+(j+1), knots_u[curU]); 01774 rowU = newRowU; 01775 } 01776 } 01777 int curV = startV+1; 01778 for(uint j=0; j<locKnotV.size()-1; j++, curV++) { 01779 if(locKnotV[j+1] != knots_v[curV]) { 01780 std::vector<double > newRowV(rowV.size()+1, 0); 01781 for(uint k=0; k<rowV.size(); k++) { 01782 #define V(x) (locKnotV[x+k]) 01783 double z = knots_v[curV]; 01784 int p = order_[1]-1; 01785 if(z < V(0) || z > V(p+1)) { 01786 newRowV[k] = rowV[k]; 01787 continue; 01788 } 01789 double alpha1 = (V(p) <= z ) ? 1 : ( z - V(0))/( V(p) - V(0)); 01790 double alpha2 = (z <= V(1)) ? 1 : ( V(p+1) - z )/( V(p+1) - V(1)); 01791 newRowV[k] += rowV[k]*alpha1; 01792 newRowV[k+1] += rowV[k]*alpha2; 01793 #undef V 01794 } 01795 locKnotV.insert(locKnotV.begin()+(j+1), knots_v[curV]); 01796 rowV= newRowV; 01797 } 01798 } 01799 std::vector<double > totalRow(fullDim, 0.0); 01800 for(uint i1=0; i1<rowU.size(); i1++) 01801 for(uint i2=0; i2<rowV.size(); i2++) 01802 totalRow[(startV+i2)*n1 + (startU+i1)] = rowV[i2]*rowU[i1]; 01803 01804 C.push_back(totalRow); 01805 } 01806 01807 if(verbose && sparseVerbose) { 01808 for(uint i=0; i<C.size(); i++) { 01809 std::cout << "|"; 01810 for(uint j=0; j<C[i].size(); j++) 01811 if(C[i][j]==0) { 01812 if(fullVerbose) 01813 std::cout << "\t"; 01814 else if(sparseVerbose) 01815 std::cout << " "; 01816 } else { 01817 if(fullVerbose) 01818 std::cout << C[i][j] << "\t"; 01819 else if(sparseVerbose) 01820 std::cout << "x"; 01821 } 01822 std::cout << "|\n"; 01823 } 01824 std::cout << std::endl; 01825 } 01826 01827 // gauss-jordan elimination 01828 int zeroColumns = 0; 01829 for(uint i=0; i<C.size() && i+zeroColumns<C[i].size(); i++) { 01830 double maxPivot = 0; 01831 int maxI = -1; 01832 for(uint j=i; j<C.size(); j++) { 01833 if(fabs(C[j][i+zeroColumns]) > maxPivot) { 01834 maxPivot = fabs(C[j][i+zeroColumns]); 01835 maxI = j; 01836 } 01837 } 01838 if(maxPivot > 0 && maxPivot < 1e-10) { 01839 C[maxI][i+zeroColumns] = 0.0; 01840 maxI = -1; 01841 } 01842 if(maxI == -1) { 01843 i--; 01844 zeroColumns++; 01845 continue; 01846 } 01847 std::vector<double> tmp = C[i]; 01848 C[i] = C[maxI]; 01849 C[maxI] = tmp; 01850 for(uint j=i+1; j<C.size(); j++) { 01851 double scale = C[j][i+zeroColumns] / C[i][i+zeroColumns]; 01852 if(scale != 0) { 01853 for(uint k=i+zeroColumns; k<C[j].size(); k++) 01854 C[j][k] -= C[i][k] * scale; 01855 } 01856 } 01857 } 01858 01859 if(verbose && sparseVerbose) { 01860 for(uint i=0; i<C.size(); i++) { 01861 std::cout << "|"; 01862 for(uint j=0; j<C[i].size(); j++) 01863 if(C[i][j]==0) { 01864 if(fullVerbose) 01865 std::cout << "\t"; 01866 else if(sparseVerbose) 01867 std::cout << " "; 01868 } else { 01869 if(fullVerbose) 01870 std::cout << C[i][j] << "\t"; 01871 else if(sparseVerbose) 01872 std::cout << "x"; 01873 } 01874 std::cout << "|\n"; 01875 } 01876 std::cout << std::endl; 01877 } 01878 01879 int rank = (nmb_bas < n1*n2-zeroColumns) ? nmb_bas : n1*n2-zeroColumns; 01880 if(verbose) { 01881 std::cout << "Matrix size : " << nmb_bas << " x " << n1*n2 << std::endl; 01882 std::cout << "Matrix rank : " << rank << std::endl; 01883 } 01884 01885 return rank == nmb_bas; 01886 } 01887 01888 double LRSplineSurface::makeIntegerKnots() { 01889 // find the smallest knot interval 01890 double smallKnotU = DBL_MAX; 01891 double smallKnotV = DBL_MAX; 01892 std::vector<double> knots_u, knots_v; 01893 getGlobalKnotVector(knots_u, knots_v); 01894 for(uint i=0; i<knots_u.size()-1; i++) 01895 if(knots_u[i+1]-knots_u[i] < smallKnotU && knots_u[i+1] != knots_u[i]) 01896 smallKnotU = knots_u[i+1]-knots_u[i]; 01897 for(uint i=0; i<knots_v.size()-1; i++) 01898 if(knots_v[i+1]-knots_v[i] < smallKnotV && knots_v[i+1] != knots_v[i]) 01899 smallKnotV = knots_v[i+1]-knots_v[i]; 01900 double scale = (smallKnotU<smallKnotV) ? smallKnotU : smallKnotV; 01901 01902 // scale all meshline values according to this 01903 Meshline *m; 01904 for(uint i=0; i<meshline_.size(); i++) { 01905 m = meshline_[i]; 01906 m->const_par_ = floor(m->const_par_/scale + 0.5); 01907 m->start_ = floor(m->start_ /scale + 0.5); 01908 m->stop_ = floor(m->stop_ /scale + 0.5); 01909 } 01910 01911 // scale all element values 01912 Element *e; 01913 for(uint i=0; i<element_.size(); i++) { 01914 e = element_[i]; 01915 e->setUmin( floor(e->umin() /scale + 0.5)); 01916 e->setVmin( floor(e->vmin() /scale + 0.5)); 01917 e->setUmax( floor(e->umax() /scale + 0.5)); 01918 e->setVmax( floor(e->vmax() /scale + 0.5)); 01919 } 01920 01921 // scale all basis functions values 01922 for(Basisfunction *b : basis_) { 01923 for(int j=0; j<order_[0]+1; j++) 01924 (*b)[0][j] = floor((*b)[0][j]/scale + 0.5); 01925 for(int j=0; j<order_[1]+1; j++) 01926 (*b)[1][j] = floor((*b)[1][j]/scale + 0.5); 01927 } 01928 01929 // scale all LRSplineSurface values 01930 start_[0] = floor(start_[0]/scale + 0.5); 01931 start_[1] = floor(start_[1]/scale + 0.5); 01932 end_[0] = floor(end_[0] /scale + 0.5); 01933 end_[1] = floor(end_[1] /scale + 0.5); 01934 01935 return scale; 01936 } 01937 01938 /************************************************************************************************************************/ 01946 std::vector<LRSplineSurface*> LRSplineSurface::getDerivativeSpace() const { 01947 int p1 = order_[0]; 01948 int p2 = order_[1]; 01949 double knotU[2*p1]; 01950 double knotV[2*p2]; 01951 for(int i=0; i<p1; i++) 01952 knotU[i] = start_[0]; 01953 for(int i=0; i<p1; i++) 01954 knotU[i+p1] = end_[0]; 01955 for(int i=0; i<p2; i++) 01956 knotV[i] = start_[1]; 01957 for(int i=0; i<p2; i++) 01958 knotV[i+p2] = end_[1]; 01959 int N1 = dim_ * (p1-1)*p2; 01960 int N2 = dim_ * p1*(p2-1); 01961 double coef1[N1]; 01962 double coef2[N2]; 01963 for(int i=0; i<N1; i++) coef1[i] = 0.0; 01964 for(int i=0; i<N2; i++) coef2[i] = 0.0; 01965 01966 LRSplineSurface *diffU = new LRSplineSurface(p1-1, p2 , p1-1, p2 , knotU+1, knotV, coef1, dim_, false); 01967 LRSplineSurface *diffV = new LRSplineSurface(p1 , p2-1, p1 , p2-1, knotU, knotV+1, coef2, dim_, false); 01968 01969 for(Meshline *m : meshline_) { 01970 diffU->insert_line(!m->span_u_line_, m->const_par_, m->start_, m->stop_, m->multiplicity_); 01971 diffV->insert_line(!m->span_u_line_, m->const_par_, m->start_, m->stop_, m->multiplicity_); 01972 } 01973 std::vector<LRSplineSurface*> results(2); 01974 results[0] = diffU; 01975 results[1] = diffV; 01976 return results; 01977 } 01978 01979 /************************************************************************************************************************/ 01987 LRSplineSurface* LRSplineSurface::getRaiseOrderSpace(int raiseOrderU, int raiseOrderV) const { 01988 int p1 = order_[0]+raiseOrderU; 01989 int p2 = order_[1]+raiseOrderV; 01990 double knotU[2*p1]; 01991 double knotV[2*p2]; 01992 for(int i=0; i<p1; i++) 01993 knotU[i] = start_[0]; 01994 for(int i=0; i<p1; i++) 01995 knotU[i+p1] = end_[0]; 01996 for(int i=0; i<p2; i++) 01997 knotV[i] = start_[1]; 01998 for(int i=0; i<p2; i++) 01999 knotV[i+p2] = end_[1]; 02000 int N = dim_ * p1*p2; 02001 double coef[N]; 02002 for(int i=0; i<N; i++) 02003 coef[i] = 0.0; 02004 LRSplineSurface *result = new LRSplineSurface(p1,p2,p1,p2,knotU, knotV, coef, dim_, false); 02005 02006 for(Meshline *m : meshline_) { 02007 int newMult = m->multiplicity_ + ((m->span_u_line_) ? raiseOrderV : raiseOrderU); 02008 result->insert_line(!m->span_u_line_, m->const_par_, m->start_, m->stop_, newMult); 02009 } 02010 02011 return result; 02012 } 02013 02014 /************************************************************************************************************************/ 02022 bool LRSplineSurface::setGlobalContinuity(int contU, int contV) { 02023 if(contU < -1 || contV < -1) 02024 return false; 02025 std::vector<Meshline*> existingLines; 02026 for(Meshline *m : meshline_) 02027 existingLines.push_back(m->copy()); 02028 02029 for(Meshline *m : existingLines) { 02030 int newMult = (m->span_u_line_) ? order_[1]-contV-1 : order_[0]-contU-1; 02031 if(newMult < 1) continue; 02032 insert_line(!m->span_u_line_, m->const_par_, m->start_, m->stop_, newMult); 02033 } 02034 02035 // clean up 02036 for(Meshline *m : existingLines) 02037 delete m; 02038 return true; 02039 } 02040 02041 /************************************************************************************************************************/ 02047 bool LRSplineSurface::decreaseContinuity(int du, int dv) { 02048 if(du < 0 || dv < 0) { 02049 return false; 02050 } 02051 std::vector<Meshline*> existingLines; 02052 for(Meshline *m : meshline_) 02053 existingLines.push_back(m->copy()); 02054 02055 for(Meshline *m : existingLines) { 02056 int newMult = m->multiplicity_ + ((m->span_u_line_) ? dv : du); 02057 int maxMult = ((m->span_u_line_) ? order_[1] : order_[0]); 02058 if(newMult > maxMult) 02059 newMult = maxMult; 02060 insert_line(!m->span_u_line_, m->const_par_, m->start_, m->stop_, newMult); 02061 } 02062 02063 // clean up 02064 for(Meshline *m : existingLines) 02065 delete m; 02066 return true; 02067 } 02068 02069 void LRSplineSurface::getSupportElements(std::vector<int> &result, const std::vector<int> &basisfunctions) const { 02070 result.clear(); 02071 std::set<int> tmp; 02072 std::vector<Element*>::iterator it; 02073 for(Basisfunction *b : basis_) { 02074 for(it=b->supportedElementBegin(); it != b->supportedElementEnd(); it++) 02075 tmp.insert((**it).getId()); 02076 } 02077 result.resize(tmp.size()); 02078 result.insert(result.begin(), tmp.begin(), tmp.end()); 02079 } 02080 02081 void LRSplineSurface::getDiagonalElements(std::vector<int> &result) const { 02082 result.clear(); 02083 for(uint i=0; i<element_.size(); i++) 02084 if(element_[i]->umin() == element_[i]->vmin()) 02085 result.push_back(i); 02086 } 02087 02088 void LRSplineSurface::getDiagonalBasisfunctions(std::vector<int> &result) const { 02089 result.clear(); 02090 int i = 0; 02091 for(Basisfunction *b : basis_) { 02092 bool isDiag = true; 02093 for(int j=0; j<order_[0]+1; j++) 02094 if((*b)[0][j] != (*b)[1][j]) 02095 isDiag = false; 02096 if(isDiag) 02097 result.push_back(i); 02098 i++; 02099 } 02100 } 02101 02102 void LRSplineSurface::getBezierElement(int iEl, std::vector<double> &controlPoints) const { 02103 controlPoints.clear(); 02104 controlPoints.resize(order_[0]*order_[1]*dim_, 0); 02105 Element *el = element_[iEl]; 02106 for(Basisfunction* b : el->support()) { 02107 std::vector<double> knotU((*b)[0].begin(), (*b)[0].end() ); 02108 std::vector<double> knotV((*b)[1].begin(), (*b)[1].end() ); 02109 int startU=-1; 02110 int startV=-1; 02111 02112 std::vector<double> rowU(1,1), rowV(1,1); 02113 02114 double min = el->umin(); 02115 double max = el->umax(); 02116 while(knotU[++startU] < min); 02117 while(true) { 02118 int p = order_[0]-1; 02119 int newI = -1; 02120 double z; 02121 if( knotU.size() < (uint) startU+order_[0] || knotU[startU+ order_[0]-1] != min) { 02122 z = min; 02123 newI = startU; 02124 } else if(knotU.size() < (uint) startU+2*order_[0] || knotU[startU+2*order_[0]-1] != max ) { 02125 z = max; 02126 newI = startU + order_[0]; 02127 } else { 02128 break; 02129 } 02130 02131 std::vector<double> newRowU(rowU.size()+1, 0); 02132 for(uint k=0; k<rowU.size(); k++) { 02133 #define U(x) ( knotU[x+k] ) 02134 if(z < U(0) || z > U(p+1)) { 02135 newRowU[k] = rowU[k]; 02136 continue; 02137 } 02138 double alpha1 = (U(p) <= z ) ? 1 : double( z - U(0)) / ( U(p) - U(0)); 02139 double alpha2 = (z <= U(1)) ? 1 : double( U(p+1) - z ) / ( U(p+1) - U(1)); 02140 newRowU[k] += rowU[k]*alpha1; 02141 newRowU[k+1] += rowU[k]*alpha2; 02142 #undef U 02143 } 02144 knotU.insert(knotU.begin()+newI, z); 02145 rowU = newRowU; 02146 } 02147 02148 min = el->vmin(); 02149 max = el->vmax(); 02150 while(knotV[++startV] < min); 02151 while(true) { 02152 int p = order_[1]-1; 02153 int newI = -1; 02154 double z; 02155 if( knotV.size() < (uint) startV+order_[1] || knotV[startV+ order_[1]-1] != min) { 02156 z = min; 02157 newI = startV; 02158 } else if(knotV.size() < (uint) startV+2*order_[1] || knotV[startV+2*order_[1]-1] != max ) { 02159 z = max; 02160 newI = startV + order_[1]; 02161 } else { 02162 break; 02163 } 02164 02165 std::vector<double> newRowV(rowV.size()+1, 0); 02166 for(uint k=0; k<rowV.size(); k++) { 02167 #define V(x) ( knotV[x+k] ) 02168 if(z < V(0) || z > V(p+1)) { 02169 newRowV[k] = rowV[k]; 02170 continue; 02171 } 02172 double alpha1 = (V(p) <= z ) ? 1 : double( z - V(0)) / ( V(p) - V(0)); 02173 double alpha2 = (z <= V(1)) ? 1 : double( V(p+1) - z ) / ( V(p+1) - V(1)); 02174 newRowV[k] += rowV[k]*alpha1; 02175 newRowV[k+1] += rowV[k]*alpha2; 02176 #undef V 02177 } 02178 knotV.insert(knotV.begin()+newI, z); 02179 rowV = newRowV; 02180 } 02181 02182 int ip = 0; 02183 for(int v=startV; v<startV+order_[1]; v++) 02184 for(int u=startU; u<startU+order_[0]; u++) 02185 for(int d=0; d<dim_; d++) 02186 controlPoints[ip++] += b->cp()[d]*rowU[u]*rowV[v]*b->w(); 02187 02188 } 02189 } 02190 02191 void LRSplineSurface::getBezierExtraction(int iEl, std::vector<double> &extractMatrix) const { 02192 Element *el = element_[iEl]; 02193 int width = order_[0]*order_[1]; 02194 int height = el->nBasisFunctions(); 02195 extractMatrix.clear(); 02196 extractMatrix.resize(width*height); 02197 02198 int rowI = 0; 02199 for(Basisfunction* b : el->support()) { 02200 int start[] = {-1,-1}; 02201 std::vector<std::vector<double> > row(2); 02202 std::vector<std::vector<double> > knot(2); 02203 for(int d=0; d<2; d++) { 02204 for(int i=0; i<order_[d]+1; i++) 02205 knot[d].push_back( (*b)[d][i] ); 02206 row[d].push_back(1); 02207 } 02208 02209 02210 for(int d=0; d<2; d++) { 02211 02212 double min = el->getParmin(d); 02213 double max = el->getParmax(d); 02214 while(knot[d][++start[d]] < min); 02215 while(true) { 02216 int p = order_[d]-1; 02217 int newI = -1; 02218 double z; 02219 if( knot[d].size() < (uint) start[d]+order_[d] || knot[d][start[d]+ order_[d]-1] != min) { 02220 z = min; 02221 newI = start[d]; 02222 } else if(knot[d].size() < (uint) start[d]+2*order_[d] || knot[d][start[d]+2*order_[d]-1] != max ) { 02223 z = max; 02224 newI = start[d] + order_[d]; 02225 } else { 02226 break; 02227 } 02228 02229 std::vector<double> newRow(row[d].size()+1, 0); 02230 for(uint k=0; k<row[d].size(); k++) { 02231 #define U(x) ( knot[d][x+k] ) 02232 if(z < U(0) || z > U(p+1)) { 02233 newRow[k] = row[d][k]; 02234 continue; 02235 } 02236 double alpha1 = (U(p) <= z ) ? 1 : double( z - U(0)) / ( U(p) - U(0)); 02237 double alpha2 = (z <= U(1)) ? 1 : double( U(p+1) - z ) / ( U(p+1) - U(1)); 02238 newRow[k] += row[d][k]*alpha1; 02239 newRow[k+1] += row[d][k]*alpha2; 02240 #undef U 02241 } 02242 knot[d].insert(knot[d].begin()+newI, z); 02243 row[d] = newRow; 02244 } 02245 02246 } 02247 02248 int colI = 0; 02249 for(int v=start[1]; v<start[1]+order_[1]; v++) 02250 for(int u=start[0]; u<start[0]+order_[0]; u++, colI++) 02251 extractMatrix[colI*height + rowI] += row[0][u]*row[1][v]*b->w(); 02252 rowI++; 02253 } 02254 } 02255 02256 void LRSplineSurface::setElementColor(double r, double g, double b) { 02257 element_red = r; 02258 element_green = g; 02259 element_blue = b; 02260 } 02261 02262 void LRSplineSurface::setBasisColor(double r, double g, double b) { 02263 basis_red = r; 02264 basis_green = g; 02265 basis_blue = b; 02266 } 02267 02268 void LRSplineSurface::setSelectedBasisColor(double r, double g, double b) { 02269 selected_basis_red = r; 02270 selected_basis_green = g; 02271 selected_basis_blue = b; 02272 } 02273 02274 void LRSplineSurface::read(std::istream &is) { 02275 start_[0] = DBL_MAX; 02276 end_[0] = -DBL_MAX; 02277 start_[1] = DBL_MAX; 02278 end_[1] = -DBL_MAX; 02279 02280 // first get rid of comments and spaces 02281 ws(is); 02282 char firstChar; 02283 char buffer[1024]; 02284 firstChar = is.peek(); 02285 while(firstChar == '#') { 02286 is.getline(buffer, 1024); 02287 ws(is); 02288 firstChar = is.peek(); 02289 } 02290 02291 // read actual parameters 02292 int nBasis, nElements, nMeshlines; 02293 is >> order_[0]; ws(is); 02294 is >> order_[1]; ws(is); 02295 is >> nBasis; ws(is); 02296 is >> nMeshlines; ws(is); 02297 is >> nElements; ws(is); 02298 is >> dim_; ws(is); 02299 is >> rational_; ws(is); 02300 02301 meshline_.resize(nMeshlines); 02302 element_.resize(nElements); 02303 std::vector<Basisfunction*> basisVector(nBasis); 02304 02305 // get rid of more comments and spaces 02306 firstChar = is.peek(); 02307 while(firstChar == '#') { 02308 is.getline(buffer, 1024); 02309 ws(is); 02310 firstChar = is.peek(); 02311 } 02312 02313 // read all basisfunctions 02314 for(int i=0; i<nBasis; i++) { 02315 Basisfunction *b = new Basisfunction(dim_, order_[0], order_[1]); 02316 b->read(is); 02317 basis_.insert(b); 02318 basisVector[i] = b; 02319 } 02320 02321 // get rid of more comments and spaces 02322 firstChar = is.peek(); 02323 while(firstChar == '#') { 02324 is.getline(buffer, 1024); 02325 ws(is); 02326 firstChar = is.peek(); 02327 } 02328 02329 for(int i=0; i<nMeshlines; i++) { 02330 meshline_[i] = new Meshline(); 02331 meshline_[i]->read(is); 02332 } 02333 02334 // get rid of more comments and spaces 02335 firstChar = is.peek(); 02336 while(firstChar == '#') { 02337 is.getline(buffer, 1024); 02338 ws(is); 02339 firstChar = is.peek(); 02340 } 02341 02342 // read elements and calculate patch boundaries 02343 for(int i=0; i<nElements; i++) { 02344 element_[i] = new Element(); 02345 element_[i]->read(is); 02346 element_[i]->updateBasisPointers(basisVector); 02347 start_[0] = (element_[i]->umin() < start_[0]) ? element_[i]->umin() : start_[0]; 02348 end_[0] = (element_[i]->umax() > end_[0] ) ? element_[i]->umax() : end_[0] ; 02349 start_[1] = (element_[i]->vmin() < start_[1]) ? element_[i]->vmin() : start_[1]; 02350 end_[1] = (element_[i]->vmax() > end_[1] ) ? element_[i]->vmax() : end_[1] ; 02351 } 02352 } 02353 02354 void LRSplineSurface::write(std::ostream &os) const { 02355 generateIDs(); 02356 os << std::setprecision(16); 02357 os << "# LRSPLINE SURFACE\n"; 02358 os << "#\tp1\tp2\tNbasis\tNline\tNel\tdim\trat\n\t"; 02359 os << order_[0] << "\t"; 02360 os << order_[1] << "\t"; 02361 os << basis_.size() << "\t"; 02362 os << meshline_.size() << "\t"; 02363 os << element_.size() << "\t"; 02364 os << dim_ << "\t"; 02365 os << rational_ << "\n"; 02366 02367 os << "# Basis functions:\n"; 02368 for(Basisfunction* b : basis_) 02369 os << *b << std::endl; 02370 os << "# Mesh lines:\n"; 02371 for(Meshline* m : meshline_) 02372 os << *m << std::endl; 02373 os << "# Elements:\n"; 02374 for(Element* e : element_) 02375 os << *e << std::endl; 02376 } 02377 02378 void LRSplineSurface::writePostscriptMesh(std::ostream &out, bool close, std::vector<int> *colorElements) const { 02379 #ifdef TIME_LRSPLINE 02380 PROFILE("Write EPS"); 02381 #endif 02382 std::vector<double> knot_u, knot_v; 02383 getGlobalUniqueKnotVector(knot_u, knot_v); 02384 double min_span_u = knot_u[1] - knot_u[0]; 02385 double min_span_v = knot_v[1] - knot_v[0]; 02386 for(uint i=1; i<knot_u.size()-1; i++) 02387 min_span_u = (min_span_u<knot_u[i+1]-knot_u[i]) ? min_span_u : knot_u[i+1]-knot_u[i]; 02388 for(uint i=1; i<knot_v.size()-1; i++) 02389 min_span_v = (min_span_v<knot_v[i+1]-knot_v[i]) ? min_span_v : knot_v[i+1]-knot_v[i]; 02390 02391 // get date 02392 time_t t = time(0); 02393 tm* lt = localtime(&t); 02394 char date[11]; 02395 sprintf(date, "%02d/%02d/%04d", lt->tm_mday, lt->tm_mon + 1, lt->tm_year+1900); 02396 02397 // get bounding box 02398 int dx = end_[0] - start_[0]; 02399 int dy = end_[1] - start_[1]; 02400 double scale = (dx>dy) ? 1000.0/dx : 1000.0/dy; 02401 // set the duplicate-knot-line (dkl) display width 02402 double dkl_range = (min_span_u>min_span_v) ? min_span_v*scale/6.0 : min_span_u*scale/6.0; 02403 int xmin = (start_[0] - dx/30.0)*scale; 02404 int ymin = (start_[1] - dy/30.0)*scale; 02405 int xmax = (end_[0] + dx/30.0)*scale + dkl_range; 02406 int ymax = (end_[1] + dy/30.0)*scale + dkl_range; 02407 02408 // print eps header 02409 out << "%!PS-Adobe-3.0 EPSF-3.0\n"; 02410 out << "%%Creator: LRSplineHelpers.cpp object\n"; 02411 out << "%%Title: LR-spline parameter domain\n"; 02412 out << "%%CreationDate: " << date << std::endl; 02413 out << "%%Origin: 0 0\n"; 02414 out << "%%BoundingBox: " << xmin << " " << ymin << " " << xmax << " " << ymax << std::endl; 02415 02416 // Fill diagonal elements when refining 02417 if(colorElements != NULL) { 02418 out << element_red << " "; 02419 out << element_green << " "; 02420 out << element_blue << " "; 02421 out << "setrgbcolor \n"; 02422 for(uint i=0; i<colorElements->size(); i++) { 02423 Element* e = element_[colorElements->at(i)]; 02424 out << "newpath\n"; 02425 out << e->umin()*scale << " " << e->vmin()*scale << " moveto\n"; 02426 out << e->umax()*scale << " " << e->vmin()*scale << " lineto\n"; 02427 out << e->umax()*scale << " " << e->vmax()*scale << " lineto\n"; 02428 out << e->umin()*scale << " " << e->vmax()*scale << " lineto\n"; 02429 out << "closepath\n"; 02430 out << "fill\n"; 02431 } 02432 } 02433 02434 out << "0 setgray\n"; 02435 out << "1 setlinewidth\n"; 02436 for(uint i=0; i<meshline_.size(); i++) { 02437 out << "newpath\n"; 02438 double dm = (meshline_[i]->multiplicity_==1) ? 0 : dkl_range/(meshline_[i]->multiplicity_-1); 02439 for(int m=0; m<meshline_[i]->multiplicity_; m++) { 02440 if(meshline_[i]->is_spanning_u()) { 02441 out << meshline_[i]->start_*scale << " " << meshline_[i]->const_par_*scale + dm*m << " moveto\n"; 02442 if(meshline_[i]->stop_ == end_[0]) 02443 out << meshline_[i]->stop_*scale+dkl_range << " " << meshline_[i]->const_par_*scale + dm*m << " lineto\n"; 02444 else 02445 out << meshline_[i]->stop_*scale << " " << meshline_[i]->const_par_*scale + dm*m << " lineto\n"; 02446 } else { 02447 out << meshline_[i]->const_par_*scale + dm*m << " " << meshline_[i]->start_*scale << " moveto\n"; 02448 if(meshline_[i]->stop_ == end_[1]) 02449 out << meshline_[i]->const_par_*scale + dm*m << " " << meshline_[i]->stop_ *scale+dkl_range << " lineto\n"; 02450 else 02451 out << meshline_[i]->const_par_*scale + dm*m << " " << meshline_[i]->stop_ *scale << " lineto\n"; 02452 } 02453 } 02454 out << "stroke\n"; 02455 } 02456 02457 if(close) 02458 out << "%%EOF\n"; 02459 } 02460 02461 void LRSplineSurface::writePostscriptElements(std::ostream &out, int nu, int nv, bool close, std::vector<int> *colorElements) const { 02462 #ifdef TIME_LRSPLINE 02463 PROFILE("Write EPS"); 02464 #endif 02465 02466 // get date 02467 time_t t = time(0); 02468 tm* lt = localtime(&t); 02469 char date[11]; 02470 sprintf(date, "%02d/%02d/%04d", lt->tm_mday, lt->tm_mon + 1, lt->tm_year+1900); 02471 02472 // get bounding box (max/min of the control points) 02473 double x[2]; 02474 double y[2]; 02475 x[0] = 1e7; 02476 x[1] = -1e7; 02477 y[0] = 1e7; 02478 y[1] = -1e7; 02479 for(Basisfunction *b : basis_) { 02480 std::vector<double>::const_iterator cp = b->cp(); 02481 x[0] = (cp[0] < x[0]) ? cp[0] : x[0]; 02482 x[1] = (cp[0] > x[1]) ? cp[0] : x[1]; 02483 y[0] = (cp[1] < y[0]) ? cp[1] : y[0]; 02484 y[1] = (cp[1] > y[1]) ? cp[1] : y[1]; 02485 } 02486 02487 double dx = x[1]-x[0]; 02488 double dy = y[1]-y[0]; 02489 double scale = (dx>dy) ? 1000.0/dx : 1000.0/dy; 02490 02491 int xmin = (x[0] - dx/20.0)*scale; 02492 int ymin = (y[0] - dy/20.0)*scale; 02493 int xmax = (x[1] + dx/20.0)*scale; 02494 int ymax = (y[1] + dy/20.0)*scale; 02495 02496 // print eps header 02497 out << "%!PS-Adobe-3.0 EPSF-3.0\n"; 02498 out << "%%Creator: LRSplineHelpers.cpp object\n"; 02499 out << "%%Title: LR-spline physical domain\n"; 02500 out << "%%CreationDate: " << date << std::endl; 02501 out << "%%Origin: 0 0\n"; 02502 out << "%%BoundingBox: " << xmin << " " << ymin << " " << xmax << " " << ymax << std::endl; 02503 02504 // Fill diagonal elements when refining 02505 if(colorElements != NULL) { 02506 out << element_red << " "; 02507 out << element_green << " "; 02508 out << element_blue << " "; 02509 out << "setrgbcolor \n"; 02510 for(uint iEl=0; iEl<element_.size(); iEl++) { 02511 bool doColor = false; 02512 for(uint j=0; j<colorElements->size(); j++) 02513 if(iEl == (uint) colorElements->at(j)) 02514 doColor = true; 02515 if(doColor) { 02516 double umin = element_[iEl]->umin(); 02517 double umax = element_[iEl]->umax(); 02518 double vmin = element_[iEl]->vmin(); 02519 double vmax = element_[iEl]->vmax(); 02520 02521 std::vector<double> pt; 02522 point(pt, umin, vmin, iEl); 02523 out << "newpath\n"; 02524 out << pt[0]*scale << " " << pt[1]*scale << " moveto\n"; 02525 for(int i=1; i<nu; i++) { // SOUTH side 02526 double u = umin + (umax-umin)*i/(nu-1); 02527 double v = vmin; 02528 point(pt, u, v, iEl, false, true); 02529 out << pt[0]*scale << " " << pt[1]*scale << " lineto\n"; 02530 } 02531 for(int i=1; i<nv; i++) { // EAST side 02532 double u = umax; 02533 double v = vmin + (vmax-vmin)*i/(nv-1); 02534 point(pt, u, v, iEl, false, false); 02535 out << pt[0]*scale << " " << pt[1]*scale << " lineto\n"; 02536 } 02537 for(int i=nu-1; i-->0; ) { // NORTH side 02538 double u = umin + (umax-umin)*i/(nu-1); 02539 double v = vmax; 02540 point(pt, u, v, iEl, true, false); 02541 out << pt[0]*scale << " " << pt[1]*scale << " lineto\n"; 02542 } 02543 for(int i=nv-1; i-->1; ) { // WEST side 02544 double u = umin; 02545 double v = vmin + (vmax-vmin)*i/(nv-1); 02546 point(pt, u, v, iEl, true, false); 02547 out << pt[0]*scale << " " << pt[1]*scale << " lineto\n"; 02548 } 02549 out << "closepath\n"; 02550 out << "fill\n"; 02551 } 02552 } 02553 } 02554 02555 out << "0 setgray\n"; 02556 out << "1 setlinewidth\n"; 02557 02558 for(uint iEl=0; iEl<element_.size(); iEl++) { 02559 double umin = element_[iEl]->umin(); 02560 double umax = element_[iEl]->umax(); 02561 double vmin = element_[iEl]->vmin(); 02562 double vmax = element_[iEl]->vmax(); 02563 02564 02565 /* DIAGONAL TEST REFINEMENT 02566 * used to color the diagonal elements (only 2x2 evaluation points) 02567 * maybe dust off this and color arbitrary elements at a later point 02568 */ 02569 /* point(corner[0], element_[iEl]->umin(), element_[iEl]->vmin(), iEl); 02570 point(corner[1], element_[iEl]->umin(), element_[iEl]->vmax(), iEl); 02571 point(corner[2], element_[iEl]->umax(), element_[iEl]->vmax(), iEl); 02572 point(corner[3], element_[iEl]->umax(), element_[iEl]->vmin(), iEl); 02573 02574 if(colorDiag && element_[iEl]->umin() == element_[iEl]->vmin()) { 02575 out << ".5 setgray\n"; 02576 if(element_[iEl]->umin() == element_[iEl]->vmin()) { 02577 out << "newpath\n"; 02578 out << corner[0][0]*scale << " " << corner[0][1]*scale << " moveto\n"; 02579 out << corner[1][0]*scale << " " << corner[1][1]*scale << " lineto\n"; 02580 out << corner[2][0]*scale << " " << corner[2][1]*scale << " lineto\n"; 02581 out << corner[3][0]*scale << " " << corner[3][1]*scale << " lineto\n"; 02582 out << "closepath\n"; 02583 out << "fill\n"; 02584 } 02585 out << "0 setgray\n"; 02586 } 02587 */ 02588 02589 std::vector<double> pt; 02590 point(pt, umin, vmin, iEl); 02591 out << "newpath\n"; 02592 out << pt[0]*scale << " " << pt[1]*scale << " moveto\n"; 02593 for(int i=1; i<nu; i++) { // SOUTH side 02594 double u = umin + (umax-umin)*i/(nu-1); 02595 double v = vmin; 02596 point(pt, u, v, iEl, false, true); 02597 out << pt[0]*scale << " " << pt[1]*scale << " lineto\n"; 02598 } 02599 for(int i=1; i<nv; i++) { // EAST side 02600 double u = umax; 02601 double v = vmin + (vmax-vmin)*i/(nv-1); 02602 point(pt, u, v, iEl, false, false); 02603 out << pt[0]*scale << " " << pt[1]*scale << " lineto\n"; 02604 } 02605 for(int i=nu-1; i-->0; ) { // NORTH side 02606 double u = umin + (umax-umin)*i/(nu-1); 02607 double v = vmax; 02608 point(pt, u, v, iEl, true, false); 02609 out << pt[0]*scale << " " << pt[1]*scale << " lineto\n"; 02610 } 02611 for(int i=nv-1; i-->1; ) { // WEST side 02612 double u = umin; 02613 double v = vmin + (vmax-vmin)*i/(nv-1); 02614 point(pt, u, v, iEl, true, false); 02615 out << pt[0]*scale << " " << pt[1]*scale << " lineto\n"; 02616 } 02617 out << "closepath\n"; 02618 out << "stroke\n"; 02619 } 02620 02621 if(close) 02622 out << "%%EOF\n"; 02623 } 02624 02625 void LRSplineSurface::writePostscriptMeshWithControlPoints(std::ostream &out, int nu, int nv) const { 02626 writePostscriptElements(out, nu, nv, false); 02627 02628 // get bounding box (max/min of the control points) 02629 double x[2]; 02630 double y[2]; 02631 x[0] = 1e7; 02632 x[1] = -1e7; 02633 y[0] = 1e7; 02634 y[1] = -1e7; 02635 for(Basisfunction *b : basis_) { 02636 std::vector<double>::const_iterator cp = b->cp(); 02637 x[0] = (cp[0] < x[0]) ? cp[0] : x[0]; 02638 x[1] = (cp[0] > x[1]) ? cp[0] : x[1]; 02639 y[0] = (cp[1] < y[0]) ? cp[1] : y[0]; 02640 y[1] = (cp[1] > y[1]) ? cp[1] : y[1]; 02641 } 02642 02643 double dx = x[1]-x[0]; 02644 double dy = y[1]-y[0]; 02645 double scale = (dx>dy) ? 1000.0/dx : 1000.0/dy; 02646 02647 double circleSize = 15.0; 02648 02649 // create the ellipse function 02650 out << "/ellipse {\n"; 02651 out << "/endangle exch def\n"; 02652 out << "/startangle exch def\n"; 02653 out << "/yrad exch def\n"; 02654 out << "/xrad exch def\n"; 02655 out << "/y exch def\n"; 02656 out << "/x exch def\n"; 02657 out << "/savematrix matrix currentmatrix def\n"; 02658 out << "x y translate\n"; 02659 out << "xrad yrad scale\n"; 02660 out << "0 0 1 startangle endangle arc\n"; 02661 out << "savematrix setmatrix\n"; 02662 out << "} def\n"; 02663 02664 // load the font to use 02665 out << "/Times-Roman findfont\n"; 02666 out << "24 scalefont\n"; 02667 out << "setfont\n"; 02668 02669 int i=-1; 02670 for(Basisfunction *b : basis_) { 02671 i++; 02672 double cp_x = b->cp(0); 02673 double cp_y = b->cp(1); 02674 // move C^{-1} text on internal functions 02675 int textX = ((*b)[0][1] == (*b)[0][order_[0]]) ? -2 : 1; 02676 int textY = ((*b)[1][1] == (*b)[1][order_[1]]) ? -2 : 1; 02677 // move text on edge functions 02678 if((*b)[0][1] == end_[0]) 02679 textX = 1; 02680 else if((*b)[0][order_[0]-1] == start_[0]) 02681 textX = -2; 02682 if((*b)[1][1] == end_[1]) 02683 textY = 1; 02684 else if((*b)[1][order_[1]-1] == start_[1]) 02685 textY = -2; 02686 02687 out << "newpath\n"; 02688 out << "0.45 0.45 0.45 setrgbcolor \n"; 02689 out << cp_x*scale << " " << cp_y*scale << " " << circleSize << " " << circleSize << " 0 360 ellipse\n"; 02690 out << "closepath fill\n"; 02691 out << "0 setgray\n"; 02692 out << cp_x*scale << " " << cp_y*scale << " " << circleSize << " " << circleSize << " 0 360 ellipse\n"; 02693 out << "closepath stroke\n"; 02694 out << "\n"; 02695 out << "newpath\n"; 02696 out << cp_x*scale + textX*circleSize << " " << cp_y*scale + textY*circleSize << " moveto\n"; 02697 out << "(" << i << ") show\n"; 02698 out << "\n"; 02699 } 02700 out << "%%EOF\n"; 02701 } 02702 02703 void LRSplineSurface::writePostscriptFunctionSpace(std::ostream &out, std::vector<int> *colorBasis, bool drawAll, bool close) const { 02704 if(drawAll) 02705 writePostscriptMesh(out, false); 02706 02707 int dx = end_[0] - start_[0]; 02708 int dy = end_[1] - start_[1]; 02709 double scale = (dx>dy) ? 1000.0/dx : 1000.0/dy; 02710 02711 double max_du = 0; 02712 double max_dv = 0; 02713 for(Basisfunction *b : basis_) { 02714 double du = (*b)[0][order_[0]] - (*b)[0][0]; 02715 double dv = (*b)[1][order_[1]] - (*b)[1][0]; 02716 max_du = (max_du > du) ? max_du : du; 02717 max_dv = (max_dv > dv) ? max_dv : dv; 02718 } 02719 02720 double scaleSize = (max_du > max_dv) ? 1.0/max_du : 1.0/max_dv; 02721 scaleSize *= 20.0; 02722 02723 // create the ellipse function 02724 out << "/ellipse {\n"; 02725 out << "/endangle exch def\n"; 02726 out << "/startangle exch def\n"; 02727 out << "/yrad exch def\n"; 02728 out << "/xrad exch def\n"; 02729 out << "/y exch def\n"; 02730 out << "/x exch def\n"; 02731 out << "/savematrix matrix currentmatrix def\n"; 02732 out << "x y translate\n"; 02733 out << "xrad yrad scale\n"; 02734 out << "0 0 1 startangle endangle arc\n"; 02735 out << "savematrix setmatrix\n"; 02736 out << "} def\n"; 02737 02738 // load the font to use 02739 out << "/Times-Roman findfont\n"; 02740 out << "24 scalefont\n"; 02741 out << "setfont\n"; 02742 02743 int i = -1; 02744 for(Basisfunction *b : basis_) { 02745 i++; 02746 double avg_u = 0; 02747 double avg_v = 0; 02748 double du = (*b)[0][order_[0]] - (*b)[0][0]; 02749 double dv = (*b)[1][order_[1]] - (*b)[1][0]; 02750 02751 // move C^{-1} text on internal functions 02752 double textOffset = 15.0; 02753 int textX = ((*b)[0][1] == (*b)[0][order_[0]]) ? -2 : 1; 02754 int textY = ((*b)[1][1] == (*b)[1][order_[1]]) ? -2 : 1; 02755 // move text on edge functions 02756 if((*b)[0][1] == end_[0]) 02757 textX = 1; 02758 else if((*b)[0][order_[0]-1] == start_[0]) 02759 textX = -2; 02760 if((*b)[1][1] == end_[1]) 02761 textY = 1; 02762 else if((*b)[1][order_[1]-1] == start_[1]) 02763 textY = -2; 02764 02765 for(int j=1; j<order_[0]; j++) 02766 avg_u += (*b)[0][j]; 02767 for(int j=1; j<order_[1]; j++) 02768 avg_v += (*b)[1][j]; 02769 avg_u /= (order_[0]-1); 02770 avg_v /= (order_[1]-1); 02771 02772 bool doColor = false; 02773 if(colorBasis != NULL) 02774 for(uint j=0; j<colorBasis->size(); j++) 02775 if(i == (uint) colorBasis->at(j)) 02776 doColor = true; 02777 02778 if(doColor) { 02779 out << selected_basis_red << " "; 02780 out << selected_basis_green << " "; 02781 out << selected_basis_blue << " "; 02782 out << "setrgbcolor \n"; 02783 } else { 02784 out << basis_red << " "; 02785 out << basis_green << " "; 02786 out << basis_blue << " "; 02787 out << "setrgbcolor \n"; 02788 } 02789 if(drawAll || doColor) { 02790 out << "newpath\n"; 02791 out << avg_u*scale << " " << avg_v*scale << " " << du*scaleSize << " " << dv*scaleSize << " 0 360 ellipse\n"; 02792 out << "closepath fill\n"; 02793 out << "0 setgray\n"; 02794 out << avg_u*scale << " " << avg_v*scale << " " << du*scaleSize << " " << dv*scaleSize << " 0 360 ellipse\n"; 02795 out << "closepath stroke\n"; 02796 } 02797 02798 out << "\n"; 02799 out << "newpath\n"; 02800 out << avg_u*scale + textX*textOffset << " " << avg_v*scale + textY*textOffset << " moveto\n"; 02801 out << "(" << i << ") show\n"; 02802 out << "\n"; 02803 } 02804 if(close) 02805 out << "%%EOF\n"; 02806 } 02807 02808 void LRSplineSurface::printElements(std::ostream &out) const { 02809 for(uint i=0; i<element_.size(); i++) { 02810 if(i<100) out << " "; 02811 if(i<10) out << " "; 02812 out << i << ": " << *element_[i] << std::endl; 02813 } 02814 } 02815 02816 #undef DOUBLE_TOL 02817 02818 } // end namespace LR 02819
1.7.6.1