/********************************************************************/ /* Copyright (c) 2017 System fugen G.K. and Yuzi Mizuno */ /* All rights reserved. */ /********************************************************************/ #include "MGCLStdAfx.h" #include "mg/Pvector.h" #include "mg/Tolerance.h" #include "mg/Box.h" #include "mg/KnotArray.h" #include "mg/SPointSeq.h" #include "mg/Straight.h" #include "mg/CCisect_list.h" #include "mg/LBRepEndC.h" #include "mg/SBRepEndC.h" #include "mg/SBRepTP.h" #include "mg/Coons.h" #include "mg/SBRep.h" #include "cskernel/Bvstan.h" using namespace std; #if defined(_DEBUG) #define new DEBUG_NEW #undef THIS_FILE static char THIS_FILE[] = __FILE__; #endif //using namespace std; //Compute 2nd derivative. MGVector get_2nd(const MGVector& first,const MGVector& second_guess,double curvature){ if(MGRZero(second_guess.len())) return second_guess; MGUnit_vector e2(second_guess); double len=first.len()*curvature; return e2*len; } //ノットをあわせた後の境界線、ブレンド関数、データポイントから点列を求める。 void bilinear_spoint_proc( const MGPvector<MGLBRep>& brepl,//ノットベクトルをあわせた後の境界線 const MGCurve& blendCrvU, //空間次元1のu方向のブレンド曲線(パラメータ、値域ともに0,1) const MGCurve& blendCrvV, //空間次元1のv方向のブレンド曲線(パラメータ、値域ともに0,1) const MGNDDArray& utau, //u方向のデータポイント const MGNDDArray& vtau, //v方向のデータポイント MGSPointSeq& spoint) //点列 { spoint = MGSPointSeq(utau.length(), vtau.length(), brepl[0]->sdim()); //4隅の点は境界線の端点の中間とする MGPosition p00 = (brepl[0]->start_point() + brepl[3]->start_point()) / 2.0, p01 = (brepl[2]->start_point() + brepl[3]->end_point()) / 2.0, p10 = (brepl[0]->end_point() + brepl[1]->start_point()) / 2.0, p11 = (brepl[1]->end_point() + brepl[2]->end_point()) / 2.0; //境界線上の点を求める int i, nu=utau.length(), nv=vtau.length(); int num1=nu-1, nvm1=nv-1; for(i=0; i<nu; i++){ spoint.store_at(i, 0, brepl[0]->eval(utau(i))); spoint.store_at(i, nvm1, brepl[2]->eval(utau(i))); } for(i=0; i<nv; i++){ spoint.store_at(0, i, brepl[3]->eval(vtau(i))); spoint.store_at(num1, i, brepl[1]->eval(vtau(i))); } //内部点を求める double spanu = brepl[0]->param_span(), spanv = brepl[1]->param_span(); for(i=1; i<num1; i++){ double blendu = (blendCrvU.eval(utau(i) / spanu))(0); double onembu=1.-blendu; MGPosition pu0(spoint(i, 0)), pu1(spoint(i, nvm1)); for(int j=1; j<nvm1; j++){ double blendv = (blendCrvV.eval(vtau(j) / spanv))(0); double onembv=1.-blendv; MGPosition p0v(spoint(0, j)), p1v(spoint(num1, j)); spoint.store_at(i,j,onembu*p0v+blendu*p1v+onembv*pu0+blendv*pu1 - (onembv*(onembu*p00+blendu*p10)+blendv*(onembu*p01+blendu*p11))); } } } //境界線、ブレンド関数、接続面を与え、対辺が同じノットベクトルのスプライン曲線 //になるように再作成して、点列、ノットベクトル、データポイントを求め、境界線の //パラメータにあわせてあった接続面をリビルドした後のパラメータに変更する。 //境界線はC1連続であり、vmin,umax,vmax,uminの順で、vmin,vmaxの向きをuminから //umaxの方向にumin,umaxの向きをvminからvmaxの方向になっているものとする。 //境界線のノットベクトルをあわせるときの誤差はline_zero()を使用している。 void bilinear_spoint( const MGCurve& blendCrvU, //空間次元1のu方向のブレンド曲線(パラメータ、値域ともに0,1) const MGCurve& blendCrvV, //空間次元1のv方向のブレンド曲線(パラメータ、値域ともに0,1) const MGPvector<MGLBRep>& perimeters, const MGSBRepTP& tp, //接続面(パラメータ範囲は境界線と同じ) MGSPointSeq& spoint, //点列 MGNDDArray& utau, //u方向のデータポイント MGNDDArray& vtau //v方向のデータポイント ){ //端末ベクトル、データポイントを求める。 MGENDCOND condk[4]={MGENDC_NO,MGENDC_NO,MGENDC_NO,MGENDC_NO}; for(int i=0; i<4; i++) if(tp.specified(i)) condk[i] = MGENDC_1D; //点列を求める MGKnotVector& knotu = perimeters[0]->knot_vector(); MGKnotVector& knotv = perimeters[1]->knot_vector(); utau = MGNDDArray(condk[3], condk[1], knotu);//std::cout<<"utau="<<utau; vtau = MGNDDArray(condk[0], condk[2], knotv);//std::cout<<"vtau="<<vtau<<endl; bilinear_spoint_proc(perimeters, blendCrvU, blendCrvV, utau, vtau, spoint); } //境界線、接続面、点列データポイント、接ベクトルから接続条件(MGSBRepEndC)を求める bool bilinear_endc( MGSPointSeq& spoint, //点列 const MGSBRepTP& tp, //接続面(パラメータ範囲は境界線と同じ) const MGPvector<MGLBRep>& perimeters, const MGCurve& blendCrvU, //空間次元1のu方向のブレンド曲線(パラメータ、値域ともに0,1) const MGCurve& blendCrvV, //空間次元1のv方向のブレンド曲線(パラメータ、値域ともに0,1) const MGNDDArray& utau, //u方向のデータポイント const MGNDDArray& vtau, //v方向のデータポイント MGSBRepEndC& endc) //端末条件 { if(utau.length() != spoint.length_u() || vtau.length() != spoint.length_v() || spoint.sdim() != 3)return false; const MGLBRep& peri0=*(perimeters[0]); const MGLBRep& peri1=*(perimeters[1]); const MGLBRep& peri2=*(perimeters[2]); const MGLBRep& peri3=*(perimeters[3]); MGENDCOND condk[4]={MGENDC_NO,MGENDC_NO,MGENDC_NO,MGENDC_NO}; const int dim=3; int m; int i,p; int lenu=utau.length(), lenv=vtau.length(); int lenum1=lenu-1, lenvm1=lenv-1; double ur[2]={utau(0),utau(lenum1)}, vr[2]={vtau(0),vtau(lenvm1)}; const int ipse[2]={1,2}; int ktp, ntp, n; const double* ttp; const double* rtp; int itanse[2]={1,1}; double tanse[6]; int irtp, ip1, ip2; double work[105], tangen[3]; int tps[4]={0,0,0,0}; for(i=0; i<4; i++) if(tp.specified(i)) tps[i]=1; if(tps[0] || tps[2]){ const MGKnotVector& knotv = peri1.knot_vector(); MGNDDArray vtau2(condk[0], condk[2], knotv);//std::cout<<"vtau="<<vtau<<endl<<"vtau2="<<vtau2<<endl; MGSPointSeq spoint2v; bilinear_spoint_proc(perimeters, blendCrvU, blendCrvV, utau, vtau2, spoint2v); spoint2v.capacity(ip1, ip2); n=vtau2.length(); double spanu = peri0.param_span(); for(m=0; m<2; m++){ //Process of perimeter num 0 and 2(v=min and max line) int perimeter=m*2; if(tps[perimeter]){ const MGLBRep& tpm = tp.TP(perimeter); ktp=tpm.order(); ntp=tpm.bdim(); ttp=tpm.knot_data(); rtp=tpm.coef_data(); irtp=tpm.line_bcoef().capacity(); //We have to generate first and 2nd derivative data for this perimeter. double v=vr[m]; MGVector deri1[2]={peri3.eval(v,1),peri1.eval(v,1)}; MGBPointSeq first(lenu,dim); for(p=0; p<dim; p++) tanse[p]=first(0,p)=deri1[0][p]; for(p=0; p<dim; p++) tanse[p+3]=first(lenum1,p)=deri1[1][p]; for(i=1; i<lenum1; i++){ double tptau=utau(i); bvstan_(ur,ktp,ntp,ttp,rtp,tptau,n,vtau2.data(), spoint2v.data(i,0,0),ipse[m],itanse,tanse,irtp, ip1,ip2,work,tangen); MGVector deri1ati(dim,tangen); first.store_at(i,deri1ati); } endc.set_1st(perimeter,first); } } } int ntemp=lenv-tps[0]-tps[2]; if((tps[0] || tps[2]) && ntemp>=2){ MGVector deri2sv[2]={peri3.eval(vr[0],2),peri1.eval(vr[0],2)}; MGVector deri2ev[2]={peri3.eval(vr[1],2),peri1.eval(vr[1],2)}; double crvtr0s=peri3.curvature(vr[0]), crvtr1s=peri1.curvature(vr[0]); double crvtr0e=peri3.curvature(vr[1]), crvtr1e=peri1.curvature(vr[1]); double spanu = peri0.param_span(); MGNDDArray taut(ntemp); taut(0)=vtau(0); for(int j=1 ; j<ntemp-1; j++) taut(j)=vtau(j+tps[0]); taut(ntemp-1)=vtau(lenv-1);//std::cout<<taut<<endl; for(i=1; i<lenum1; i++){ double tptau=utau(i); MGLBRepEndC endc0, endc2; MGBPointSeq bpt(ntemp,3); bpt.store_at(0,spoint(i,0)); double blendu = (blendCrvU.eval(tptau / spanu))(0); if(tps[0]){ MGVector deri1=MGVector(3,endc.first(0)(i)); endc0.set_1st(deri1); MGVector deri2ati0=deri2sv[0].interpolate_by_rotate(blendu,deri2sv[1]); deri2ati0=get_2nd(deri1,deri2ati0,crvtr0s+(crvtr1s-crvtr0s)*blendu); endc0.set_2nd(deri2ati0); } for(int j=1 ; j<ntemp-1; j++) bpt.store_at(j,spoint(i,j+tps[0])); if(tps[2]){ MGVector deri1=MGVector(3,endc.first(2)(i)); endc2.set_1st(deri1); MGVector deri2ati0=deri2ev[0].interpolate_by_rotate(blendu,deri2ev[1]); deri2ati0=get_2nd(deri1,deri2ati0,crvtr0e+(crvtr1e-crvtr0e)*blendu); endc2.set_2nd(deri2ati0); } bpt.store_at(ntemp-1,spoint(i,lenv-1)); int k=6; int new_bdim=ntemp+tps[0]*2+tps[2]*2; if(new_bdim<k) k=new_bdim; MGLBRep lbt(k,endc0,endc2,taut,bpt); if(tps[0]) spoint.store_at(i,1,lbt.eval(vtau(1))); if(tps[2]) spoint.store_at(i,lenv-2,lbt.eval(vtau(lenv-2))); } } if(tps[3] || tps[1]){ const MGKnotVector& knotu = peri0.knot_vector(); MGNDDArray utau2(condk[3], condk[1], knotu);//std::cout<<"utau="<<utau<<endl<<"utau2="<<utau2<<endl; MGSPointSeq spoint2u; bilinear_spoint_proc(perimeters, blendCrvU, blendCrvV, utau2, vtau, spoint2u); int psizeu, psizev; spoint2u.capacity(psizeu, psizev); n=utau2.length(); ip1=1; ip2=psizeu*psizev; const int iv[2]={3,1}; double spanv = peri1.param_span(); for(m=0; m<2; m++){ //Process of perimeter num 3 and 1(u=min and max line) int perimeter=iv[m]; if(tps[perimeter]){ const MGLBRep& tpm = tp.TP(perimeter); ktp=tpm.order(); ntp=tpm.bdim(); ttp=tpm.knot_data(); rtp=tpm.coef_data(); irtp=tpm.line_bcoef().capacity(); //We have to generate first and 2nd derivative data for this perimeter. double u=ur[m]; MGVector deri1[2]={peri0.eval(u,1),peri2.eval(u,1)}; MGBPointSeq first(lenv,dim); for(p=0; p<dim; p++) tanse[p]=first(0,p)=deri1[0][p]; for(p=0; p<dim; p++) tanse[p+3]=first(lenvm1,p)=deri1[1][p]; for(i=1; i<lenvm1; i++){ double tptau=vtau(i); bvstan_(vr,ktp,ntp,ttp,rtp,tptau,n,utau2.data(), spoint2u.data(0,i,0),ipse[m],itanse,tanse,irtp, ip1,ip2,work,tangen); MGVector deri1ati(dim,tangen); first.store_at(i,deri1ati); } endc.set_1st(perimeter,first); } } } ntemp=lenu-tps[3]-tps[1]; if((tps[3] || tps[1]) && ntemp>=2){ MGVector deri2su[2]={peri0.eval(ur[0],2),peri2.eval(ur[0],2)}; MGVector deri2eu[2]={peri0.eval(ur[1],2),peri2.eval(ur[1],2)}; double crvtr0s=peri0.curvature(ur[0]), crvtr1s=peri2.curvature(ur[0]); double crvtr0e=peri0.curvature(ur[1]), crvtr1e=peri2.curvature(ur[1]); double spanv = peri1.param_span(); MGNDDArray taut(ntemp); taut(0)=utau(0); for(int j=1 ; j<ntemp-1; j++) taut(j)=utau(j+tps[3]); taut(ntemp-1)=utau(lenu-1);//std::cout<<taut<<endl; for(i=1; i<lenvm1; i++){ double tptau=vtau(i); MGLBRepEndC endc0, endc2; MGBPointSeq bpt(ntemp,3); bpt.store_at(0,spoint(0,i)); double blendv = (blendCrvV.eval(tptau / spanv))(0); if(tps[3]){ MGVector deri2ati0=deri2su[0].interpolate_by_rotate(blendv,deri2su[1]); endc0.set_1st(MGVector(3,endc.first(3)(i))); endc0.set_2nd(deri2ati0); } for(int j=1 ; j<ntemp-1; j++) bpt.store_at(j,spoint(j+tps[3],i)); if(tps[1]){ MGVector deri2ati0=deri2eu[0].interpolate_by_rotate(blendv,deri2eu[1]); endc2.set_1st(MGVector(3,endc.first(1)(i))); endc2.set_2nd(deri2ati0); } bpt.store_at(ntemp-1,spoint(lenu-1,i)); int k=6; int new_bdim=ntemp+tps[3]*2+tps[1]*2; if(new_bdim<k) k=new_bdim; MGLBRep lbt(k,endc0,endc2,taut,bpt); if(tps[3]) spoint.store_at(1,i,lbt.eval(utau(1))); if(tps[1]) spoint.store_at(lenu-2,i,lbt.eval(utau(lenu-2))); } } return true; } double get_deri_coef(double t0, double t1,double alpha, double t){ if(t<=.5) return (2.*t*(1.-alpha)+2.*alpha*t0-(t0+t1))/(t0-t1); return (2.*t*(1.-alpha)+2.*alpha*t1-(t0+t1))/(t1-t0); } void get_derivatives( bool along_u, //true if for perimeter 0 and 2. //false if for perimeter 3 and 1. const MGSBRep& surf, //Original surface. const MGSBRepTP& tp, //接続面(パラメータ範囲は境界線と同じ) MGPvector<MGLBRep>& derivatives, //array of derivatives[4], must have size 4. const double* alpha //derivative magnitude coefficients for //perimeter i in alpha[i]. //=1. is normal, if<1., curvature will be made large, and //if>1. curvature will be made small. ){ int sd=surf.sdim(); const MGKnotVector& tu=surf.knot_vector_u(); const MGKnotVector& tv=surf.knot_vector_v(); if(along_u){ int nu=surf.bdim_u(); MGNDDArray utau; utau.update_from_knot(tu); MGBPointSeq deris_u(nu,sd); double v[2]={tv.param_s(),tv.param_e()}; double u0=tu.param_s(), u1=tu.param_e(); for(int m=0; m<2; m++){//for perimeter 0 and 2. int perim=m*2; double vm=v[m]; for(int i=0; i<nu; i++){ double utaui=utau[i]; MGVector deri=surf.eval(utaui,vm,0,1);//std::cout<<deri; if(tp.specified(perim)){ MGVector N=tp.TP(perim).eval(utaui);//Normal of tangent plane. double vlen=deri.len(); deri=MGUnit_vector(N*deri*N)*vlen;//std::cout<<deri<<endl; } if(alpha) deri*=get_deri_coef(u0,u1,alpha[perim],utaui); deris_u.store_at(i,deri); } derivatives.reset(perim,new MGLBRep(utau,deris_u,tu)); //std::cout<<"driv at perim="<<perim<<","<<*(derivatives[perim])<<endl; } }else{ int nv=surf.bdim_v(); MGNDDArray vtau; vtau.update_from_knot(tv); MGBPointSeq deris_v(nv,sd); double u[2]={tu.param_s(),tu.param_e()}; double v0=tv.param_s(), v1=tv.param_e(); int peri[2]={3,1}; for(int m=0; m<2; m++){//for perimeter 3 and 1. int perim=peri[m]; double um=u[m]; for(int j=0; j<nv; j++){ double vtauj=vtau[j]; MGVector deri=surf.eval(um,vtauj,1,0);//std::cout<<deri; if(tp.specified(perim)){ MGVector N=tp.TP(perim).eval(vtauj); //Normal of tangent plane at surf(um,vtau[j]). double vlen=deri.len(); deri=MGUnit_vector(N*deri*N)*vlen;//std::cout<<deri<<endl; } if(alpha) deri*=get_deri_coef(v0,v1,alpha[perim], vtauj); deris_v.store_at(j,deri); } derivatives.reset(perim,new MGLBRep(vtau,deris_v,tv)); //std::cout<<"driv at perim="<<perim<<","<<*(derivatives[perim])<<endl; } } } //境界線、ブレンド関数、接続面を与え、対辺が同じノットベクトルのスプライン曲線 //になるように再作成して、点列、ノットベクトル、データポイントを求め、境界線の //パラメータにあわせてあった接続面をリビルドした後のパラメータに変更する。 //境界線はC1連続であり、vmin,umax,vmax,uminの順で、vmin,vmaxの向きをuminから //umaxの方向にumin,umaxの向きをvminからvmaxの方向になっているものとする。 //境界線のノットベクトルをあわせるときの誤差はline_zero()を使用している。 void rebuild_curve( const MGCurve* edge_crvl[4], //境界線リスト(vmin,umax,vmax,uminの順) MGSBRepTP& tp, //接続面(パラメータ範囲は境界線と同じ) MGPvector<MGLBRep>& perimeters ){ //境界線がC1連続か、接続面のパラメータ範囲が境界線と同じかどうか調べる for(int i=0; i<4; i++){ //if(!edge_crvl[i]->cn_continuity(1)) return -1; //C1連続性のチェック if(!tp.specified(i)) continue; MGLBRep& tpi=tp.TP(i); const MGInterval cirange=edge_crvl[i]->param_range(); if(cirange != tpi.param_range()) tpi.change_range(cirange.low_point(),cirange.high_point()); } //Make the tolerances smaller. double lzero=MGTolerance::line_zero();//save line zero; MGTolerance::set_line_zero(lzero*.3); double azero=MGTolerance::angle_zero();//save line zero; MGTolerance::set_angle_zero(azero*.3); //リビルドを行う int k=4; //オーダーは4とする std::vector<const MGCurve*> temp_crvl(2); perimeters.resize(4); temp_crvl[0] = edge_crvl[0]; temp_crvl[1] = edge_crvl[2]; MGLBRep* tplb[2]={0,0}; if(tp.specified(0)) tplb[0]=&(tp.TP(0)); if(tp.specified(2)) tplb[1]=&(tp.TP(2)); MGPvector<MGLBRep> temp_brepl = rebuild_knot(temp_crvl, k, tplb); //if(!temp_brepl.size()){ // MGTolerance::set_line_zero(lzero);return -3; //} tp.set_TP(0,std::unique_ptr<MGLBRep>(tplb[0])); tp.set_TP(2,std::unique_ptr<MGLBRep>(tplb[1])); perimeters.reset(0, temp_brepl.release(0)); perimeters.reset(2, temp_brepl.release(1)); temp_crvl[0] = edge_crvl[1]; temp_crvl[1] = edge_crvl[3]; tplb[0]=tplb[1]=0; if(tp.specified(1)) tplb[0]=&(tp.TP(1)); if(tp.specified(3)) tplb[1]=&(tp.TP(3)); MGPvector<MGLBRep> temp_brepl2 = rebuild_knot(temp_crvl, k, tplb); //if(!temp_brepl2.size()){ // MGTolerance::set_line_zero(lzero);return -3; //} tp.set_TP(1,std::unique_ptr<MGLBRep>(tplb[0])); tp.set_TP(3,std::unique_ptr<MGLBRep>(tplb[1])); perimeters.reset(1, temp_brepl2.release(0)); perimeters.reset(3, temp_brepl2.release(1)); //Restore the original tolerance. MGTolerance::set_line_zero(lzero); MGTolerance::set_angle_zero(azero); MGLBRep& l0=*(perimeters[0]); MGLBRep& l1=*(perimeters[1]); MGLBRep& l2=*(perimeters[2]); MGLBRep& l3=*(perimeters[3]); //make the knots fine. MGKnotVector tu(l0.knot_vector()), tv(l1.knot_vector()); tu.change_knot_number(tu.bdim()*3); tv.change_knot_number(tv.bdim()*3); MGCurve* crv2[4]; for(int i=0; i<4; i++) crv2[i]=edge_crvl[i]->clone(); double u0=tu.param_s(), u1=tu.param_e(); double v0=tv.param_s(), v1=tv.param_e(); crv2[0]->change_range(u0,u1); crv2[1]->change_range(v0,v1); crv2[2]->change_range(u0,u1); crv2[3]->change_range(v0,v1); l0=MGLBRep(*(crv2[0]),tu); l1=MGLBRep(*(crv2[1]),tv); l2=MGLBRep(*(crv2[2]),tu); l3=MGLBRep(*(crv2[3]),tv); for(int i=0; i<4; i++){ delete crv2[i]; if(!tp.specified(i)) continue; MGKnotVector* t=&tu; if(i%2) t=&tv; MGLBRep& tpi=tp.TP(i); tpi=MGLBRep(tpi,*t); } //Adjust the corner points. MGPosition P0=(l0.start_point()+l3.start_point())*.5; MGPosition P1=(l0.end_point()+l1.start_point())*.5; MGPosition P2=(l1.end_point()+l2.end_point())*.5; MGPosition P3=(l3.end_point()+l2.start_point())*.5; double fixp[2]; //P0 fixp[0]=l0.param_e(); l0.move(2,l0.param_s(),P0,fixp); fixp[0]=l3.param_e(); l3.move(2,l3.param_s(),P0,fixp); //P1 fixp[0]=l0.param_s(); l0.move(2,l0.param_e(),P1,fixp); fixp[0]=l1.param_e(); l1.move(2,l1.param_s(),P1,fixp); //P2 fixp[0]=l1.param_s(); l1.move(2,l1.param_e(),P2,fixp); fixp[0]=l2.param_s(); l2.move(2,l2.param_e(),P2,fixp); //P3 fixp[0]=l2.param_e(); l2.move(2,l2.param_s(),P3,fixp); fixp[0]=l3.param_s(); l3.move(2,l3.param_e(),P3,fixp); } void rebuild_curve( MGPvector<MGLBRep>& edges, MGSBRepTP& tp, //接続面(パラメータ範囲は境界線と同じ) MGPvector<MGLBRep>& perimeters ){ const MGCurve* edge_crvl[4]; for(int i=0; i<4; i++) edge_crvl[i]=edges[i]; rebuild_curve(edge_crvl,tp,perimeters); } //4本の境界線、ブレンド関数、接続面を与えて面を生成する。 //境界線はvmin,umax,vmax,uminの順で、vmin,vmaxの向きをuminからumaxの方向に //umin,umaxの向きをvminからvmaxの方向になっているものとする。境界線のノットベクトル //をあわせるときの誤差はline_zero()を使用している。ブレンド曲線はパラメータ範囲0,1 //で値域も0,1である。接続面(MGSBRepTP)のパラメータ範囲は各境界線と同じとする。 // エラーコード: // 0: 正常終了 // -1: 境界線がC1連続でなかった // -2: 接続面のパラメータ範囲が境界線と違った // -3: 境界線のノットベクトルをあわせられなかった // -4: 接続面のパラメータ範囲を変更できなかった // -5: 点列が求まらなかった // -6: 端末条件が求まらなかった // -7: 面が生成できなかった MGSBRep::MGSBRep( const MGCurve* edge_crvl[4], //境界線リスト(vmin,umax,vmax,uminの順,辺番号0,1,2,3の順) const MGCurve& blendCrvU, //空間次元1のu方向のブレンド曲線(パラメータ、値域ともに0,1) const MGCurve& blendCrvV, //空間次元1のv方向のブレンド曲線(パラメータ、値域ともに0,1) const MGSBRepTP& tp, //接続面(パラメータ範囲は境界線と同じ) int& error):MGSurface() //エラーコード { //点列、ノットベクトル、データポイントを求める //for(int iii=0;iii<4; iii++) std::cout<<(*edge_crvl[iii]); MGSBRepTP tempTP(tp); MGSPointSeq spoint; MGNDDArray utau, vtau; MGPvector<MGLBRep> perimeters; //1. Adjust parameter range of the tp and edge_crvl. rebuild_curve(edge_crvl,tempTP,perimeters); //2. Build spoint data. bilinear_spoint(blendCrvU,blendCrvV,perimeters,tempTP,spoint,utau,vtau); //3. 接続条件(MGSBRepEndC)を求める MGSBRepEndC endc; if(!bilinear_endc(spoint, tempTP,perimeters, blendCrvU, blendCrvV, utau, vtau, endc)){ error = -6; return; } //4. 点列、データポイント、ノットベクトル、接続条件より面を生成する MGKnotVector& knotu=perimeters[0]->knot_vector(); MGKnotVector& knotv=perimeters[1]->knot_vector(); //std::cout<<utau<<knotu<<vtau<<knotv<<endl;////////////////////////////////// *this = MGSBRep(endc, utau, vtau, spoint, knotu, knotv); //if(error)error = -7; } //Easy to use version of the above constructor. //When blendCCrvU,V were null, straight line from 0. to 1. will be used. MGSBRep::MGSBRep( MGPvector<MGLBRep>& edges, //境界線リスト(vmin,umax,vmax,uminの順,辺番号0,1,2,3の順) int& error , //エラーコードが出力される const MGCurve* blendCrvU,//空間次元1のu方向のブレンド曲線(パラメータ、値域ともに0,1) const MGCurve* blendCrvV //空間次元1のv方向のブレンド曲線(パラメータ、値域ともに0,1) ){ MGPosition zero(1); zero(0)=0.; MGPosition one(1); one(0)=1.; MGStraight zero_one(one,zero); const MGCurve* blendu; const MGCurve* blendv; if(blendCrvU) blendu=blendCrvU; else blendu=&zero_one; if(blendCrvV) blendv=blendCrvV; else blendv=&zero_one; MGSBRepTP tp; MGSPointSeq spoint; MGNDDArray utau, vtau; MGPvector<MGLBRep> perimeters; //1. Adjust parameter range of the edge_crvl. rebuild_curve(edges,tp,perimeters); //2. Build spoint data. bilinear_spoint(*blendu,*blendv,perimeters,tp,spoint,utau,vtau); //3. 点列、データポイント、ノットベクトル、接続条件より面を生成する MGKnotVector& knotu=perimeters[0]->knot_vector(); MGKnotVector& knotv=perimeters[1]->knot_vector(); //std::cout<<utau<<knotu<<vtau<<knotv<<endl;////////////////////////////////// MGSBRepEndC endc; *this = MGSBRep(endc, utau, vtau, spoint, knotu, knotv); } MGSBRep::MGSBRep( const MGCurve* edge_crvl[4], //境界線リスト(vmin,umax,vmax,uminの順,辺番号0,1,2,3の順) const MGSBRepTP& tp, //接続面(パラメータ範囲は境界線と同じ) int& error, //エラーコード const double* alpha //derivative magnitude coefficients for perimeter i in alpha[i]. //=1. is normal, if<1., curvature will be made large, and //if>1. curvature will be made small. //If alpha=null, all of the coefs are assumed to be 1. ):MGSurface(){ // std::cout<<*(edge_crvl[0])<<endl; //std::cout<<MGTolerance::instance(); MGSBRepTP tp2(tp); MGPvector<MGLBRep> perimeters; rebuild_curve(edge_crvl,tp2,perimeters); //std::cout<<tp2<<endl; //Save the original parameter range. MGKnotVector& tu=perimeters[0]->knot_vector(); double u0=tu.param_s(), u1=tu.param_e(); MGKnotVector& tv=perimeters[1]->knot_vector(); double v0=tv.param_s(), v1=tv.param_e(); //Change parameter range to (0., 1.). int i; for(i=0; i<4; i++){ perimeters[i]->change_range(0., 1.); if(tp2.specified(i)) tp2.TP(i).change_range(0., 1.); } //std::cout<<tp2<<endl; MGPvector<MGLBRep> derivatives(4); bool along_u=true; //Method 1. MGSBRep ruled0(along_u,perimeters);//std::cout<<ruled0; MGSBRep ruled1(!along_u,perimeters);//std::cout<<ruled1; // get_derivatives(along_u,ruled0,tp2,derivatives,alpha); // get_derivatives(!along_u,ruled1,tp2,derivatives,alpha); //Method 2. MGSPointSeq sp(ruled0.surface_bcoef()+ruled1.surface_bcoef()); sp*=.5; MGSBRep ruled01(sp,ruled0.knot_vector_u(), ruled0.knot_vector_v());//std::cout<<ruled01; get_derivatives(along_u,ruled01,tp2,derivatives,alpha); get_derivatives(!along_u,ruled01,tp2,derivatives,alpha); //Method 3. /* MGSBRep bilinear(perimeters,error); bilinear.change_range(1,0.,1.); bilinear.change_range(0,0.,1.); get_derivatives(along_u,bilinear,tp2,derivatives,alpha); get_derivatives(!along_u,bilinear,tp2,derivatives,alpha); */ double taumax[4]; /* double cosmax[4]; tp2.get_perimeters_max_cos(derivatives,taumax,cosmax); for(i=0; i<4; i++) std::cout<<endl<<"i="<<i<<":taumax="<<taumax[i]<<", cosmax="<<cosmax[i];std::cout<<endl; */ //Save the old knot vector. MGKnotVector tu2(tu), tv2(tv); //Construct a coons' patch. MGCoons coons(perimeters,derivatives); MGNDDArray utau(3,tu2.bdim()-2,tu2); tv2.change_knot_number(tv2.bdim()*3); MGNDDArray vtau(3,tv2.bdim()-2,tv2); MGSPointSeq spoint(utau.length(),vtau.length(),3); coons.eval(utau,vtau,spoint); MGSBRepEndC endc(utau,vtau,coons); MGSBRep* surf=new MGSBRep(endc,utau,vtau,spoint,tu2,tv2); //if(error){ error = -7; return ;} double sinmax[4]; double error_max=tp2.get_perimeters_max_sin(*surf,taumax,sinmax); error_max*=1.005; double azero=MGTolerance::angle_zero(); if(error_max<azero) error_max=azero; double line0=MGTolerance::line_zero(); double line02=line0*2.; MGSBRep* surf2=0; //Remove knot by the minimum line zero that satisfy the angle continuity. for(int q=0; q<4; q++){ delete surf2; surf2=new MGSBRep(*surf); line02*=.5; MGTolerance::set_line_zero(line02); surf2->remove_knot(); double max2=tp2.get_perimeters_max_sin(*surf2,taumax,sinmax); if(max2<=error_max) break; } delete surf; surf=surf2; MGTolerance::set_line_zero(line0); surf->change_range(1,u0,u1); surf->change_range(0,v0,v1); *this=*surf; delete surf; }