|
LR-splines
0.5
|
00001 #include <stdio.h> 00002 #include <iostream> 00003 #include <string.h> 00004 #include <fstream> 00005 #include <sstream> 00006 #ifdef HAS_GOTOOLS 00007 #include <GoTools/geometry/SplineSurface.h> 00008 #include <GoTools/trivariate/SplineVolume.h> 00009 #include <GoTools/geometry/ObjectHeader.h> 00010 #endif 00011 #include "LRSpline/LRSplineSurface.h" 00012 #include "LRSpline/LRSplineVolume.h" 00013 #include "LRSpline/Profiler.h" 00014 #include "LRSpline/Element.h" 00015 #include "LRSpline/MeshRectangle.h" 00016 00017 using namespace LR; 00018 using namespace std; 00019 00020 int main(int argc, char **argv) { 00021 #ifdef TIME_LRSPLINE 00022 Profiler prof(argv[0]); 00023 #endif 00024 00025 // set default parameter values 00026 int p1 = 3; 00027 int p2 = 2; 00028 int p3 = 4; 00029 int n1 = 6; 00030 int n2 = 5; 00031 int n3 = 8; 00032 int it = 9; 00033 int dim = 3; 00034 bool rat = false; 00035 bool vol = false; 00036 char *lrInitMesh = NULL; 00037 stringstream parameters; 00038 parameters << " parameters: \n" \ 00039 " -p1 <n> polynomial ORDER (degree+1) in first parametric direction\n" \ 00040 " -p2 <n> polynomial order in second parametric direction\n" \ 00041 " -p3 <n> polynomial order in third parametric direction (only trivariate volumes)\n" \ 00042 " -p <n> polynomial order in all parametric directions \n" \ 00043 " -n1 <n> number of basis functions in first parametric direction\n" \ 00044 " -n2 <n> number of basis functions in second parametric direction\n" \ 00045 " -n3 <n> number of basis functions in third parametric direction (only trivariate volumes)\n" \ 00046 " -n <n> number of basis functions in all parametric directions\n" \ 00047 " -it <n> number of evaluation points per element\n" \ 00048 " -in: <s> make the LRSplineSurface <s> the initial mesh\n"\ 00049 " -help display (this) help screen\n"; 00050 parameters << " default values\n"; 00051 parameters << " -p = { " << p1 << ", " << p2 << ", " << p3 << " }\n"; 00052 parameters << " -n = { " << n1 << ", " << n2 << ", " << n3 << " }\n"; 00053 parameters << " -it = { " << it << " }\n"; 00054 parameters << " -vol = { " << ((vol)?"true":"false") << " }\n"; 00055 00056 // read input 00057 for(int i=1; i<argc; i++) { 00058 if(strcmp(argv[i], "-p1") == 0) { 00059 p1 = atoi(argv[++i]); 00060 } else if(strcmp(argv[i], "-p2") == 0) { 00061 p2 = atoi(argv[++i]); 00062 } else if(strcmp(argv[i], "-p3") == 0) { 00063 p3 = atoi(argv[++i]); 00064 vol = true; 00065 } else if(strcmp(argv[i], "-p") == 0) { 00066 p1 = p2 = p3 = atoi(argv[++i]); 00067 } else if(strcmp(argv[i], "-n1") == 0) { 00068 n1 = atoi(argv[++i]); 00069 } else if(strcmp(argv[i], "-n2") == 0) { 00070 n2 = atoi(argv[++i]); 00071 } else if(strcmp(argv[i], "-n3") == 0) { 00072 n3 = atoi(argv[++i]); 00073 vol = true; 00074 } else if(strcmp(argv[i], "-n") == 0) { 00075 n1 = n2 = n3 = atoi(argv[++i]); 00076 } else if(strcmp(argv[i], "-it") == 0) { 00077 it = atoi(argv[++i]); 00078 } else if(strncmp(argv[i], "-in:", 4) == 0) { 00079 lrInitMesh = argv[i]+4; 00080 } else if(strcmp(argv[i], "-vol") == 0) { 00081 vol = true; 00082 } else if(strcmp(argv[i], "-help") == 0) { 00083 cout << "usage: " << argv[0] << "[parameters] " << endl << parameters.str(); 00084 exit(0); 00085 } else { 00086 cerr << "usage: " << argv[0] << "[parameters] " << endl << parameters.str(); 00087 exit(1); 00088 } 00089 } 00090 00091 // do some error testing on input 00092 if(n1 < p1) { 00093 cerr << "ERROR: n1 must be greater or equal to p1\n"; 00094 exit(2); 00095 } else if(n2 < p2) { 00096 cerr << "ERROR: n2 must be greater or equal to p2\n"; 00097 exit(2); 00098 } else if(n3 < p3) { 00099 cerr << "ERROR: n3 must be greater or equal to p3\n"; 00100 exit(2); 00101 } 00102 00103 LRSplineSurface *lrs; 00104 LRSplineVolume *lrv; 00105 LRSpline *lr; 00106 00107 if(lrInitMesh == NULL) { 00108 // make a uniform integer knot vector 00109 double knot_u[n1+p1]; 00110 double knot_v[n2+p2]; 00111 double knot_w[n3+p3]; 00112 for(int i=0; i<p1+n1; i++) 00113 knot_u[i] = (i<p1) ? 0 : (i>n1) ? n1-p1+1 : i-p1+1; 00114 for(int i=0; i<p2+n2; i++) 00115 knot_v[i] = (i<p2) ? 0 : (i>n2) ? n2-p2+1 : i-p2+1; 00116 for(int i=0; i<p3+n3; i++) 00117 knot_w[i] = (i<p3) ? 0 : (i>n3) ? n3-p3+1 : i-p3+1; 00118 00119 // create a list of random control points (all between 0.1 and 1.1) 00120 int nCP = (vol) ? n1*n2*n3 : n1*n2; 00121 nCP *= (dim+rat); 00122 double cp[nCP]; 00123 int k=0; 00124 for(int i=0; i<nCP; i++) // 839 as a generator over Z_853 gives a period of 425. Should suffice 00125 cp[k++] = (i*839 % 853) / 853.0 + 0.1; // rational weights also random and thus we need >0 00126 00127 // make two spline objects 00128 if(vol) { 00129 lrv = new LRSplineVolume(n1, n2, n3, p1, p2, p3, knot_u, knot_v, knot_w, cp, dim, rat); 00130 lr = lrv; 00131 } else { 00132 lrs = new LRSplineSurface(n1, n2, p1, p2, knot_u, knot_v, cp, dim, rat); 00133 lr = lrs; 00134 } 00135 } else { 00136 #ifdef HAS_GOTOOLS 00137 ifstream inputfile; 00138 inputfile.open(lrInitMesh); 00139 if(!inputfile.is_open()) { 00140 cerr << "Error: could not open file " << lrInitMesh << endl; 00141 exit(3); 00142 } 00143 Go::ObjectHeader head; 00144 Go::SplineSurface *ss = new Go::SplineSurface(); 00145 Go::SplineVolume *sv = new Go::SplineVolume(); 00146 inputfile >> head; 00147 if(head.classType() == Go::Class_SplineVolume) { 00148 vol = true; 00149 inputfile >> *sv; 00150 lr = lrv = new LRSplineVolume(sv); 00151 } else if(head.classType() == Go::Class_SplineSurface) { 00152 vol = false; 00153 inputfile >> *ss; 00154 lr = lrs = new LRSplineSurface(ss); 00155 } 00156 #endif 00157 } 00158 lr->generateIDs(); 00159 00160 00161 // ---------------- Do evaluation on known elements -------------- 00162 vector<vector<double> > result; 00163 vector<double> pt; 00164 int maxDerivs = max(p1,p2); 00165 if(vol) maxDerivs = max(maxDerivs, p3); 00166 maxDerivs -= 2; 00167 char string[256]; 00168 bool firstPrint = true; 00169 for(Element *el : lr->getAllElements()) { 00170 double u = (el->getParmin(0) + el->getParmax(0))/2.0; 00171 double v = (el->getParmin(1) + el->getParmax(1))/2.0; 00172 double w = (vol) ? (el->getParmin(2) + el->getParmax(2))/2.0 : 0.0; 00173 for(int i=0; i<it; i++) { 00174 { 00175 PROFILE("w/id basis"); 00176 if(vol) 00177 lrv->computeBasis(u,v,w, result, 0, el->getId()); 00178 else 00179 lrs->computeBasis(u,v, result, 0, el->getId()); 00180 } 00181 { 00182 sprintf(string, "w/id %d derivs", maxDerivs); 00183 PROFILE(string); 00184 if(vol) 00185 lrv->computeBasis(u,v,w, result, maxDerivs, el->getId()); 00186 else 00187 lrs->computeBasis(u,v, result, maxDerivs, el->getId()); 00188 } 00189 if(firstPrint && i==0) 00190 cout << "w/id derivs result size: " << result.size() << " x " << result[0].size() << endl; 00191 { 00192 PROFILE("w/id point"); 00193 if(vol) 00194 lrv->point(pt, u,v,w, el->getId()); 00195 else 00196 lrs->point(pt, u,v, el->getId()); 00197 } 00198 { 00199 PROFILE("wo/id basis"); 00200 if(vol) 00201 lrv->computeBasis(u,v,w, result, 0); 00202 else 00203 lrs->computeBasis(u,v, result, 0); 00204 } 00205 { 00206 sprintf(string, "wo/id %d derivs", maxDerivs); 00207 PROFILE(string); 00208 if(vol) 00209 lrv->computeBasis(u,v,w, result, maxDerivs); 00210 else 00211 lrs->computeBasis(u,v, result, maxDerivs); 00212 } 00213 if(firstPrint && i==0) 00214 cout << "wo/id derivs result size: " << result.size() << " x " << result[0].size() << endl; 00215 { 00216 PROFILE("wo/id point"); 00217 if(vol) 00218 lrv->point(pt, u,v,w); 00219 else 00220 lrs->point(pt, u,v ); 00221 } 00222 } 00223 firstPrint = false; 00224 } 00225 00226 }
1.7.6.1