1 /** @file gsThbs_geometry_test.cpp
2
3 @brief Tests for THB.
4
5 Random tests with:
6 - random hierarchical subdomain configuration (generated by randomized box insertions to the Quadtree)
7 - uniform knot vectors with different multiplicities- multiplicity added by random. max multiplicity is equal to degree
8
9 Parameter to change for different test cases:
10 T_nlevels - the maximum numbre of levels
11 test_num - number of tests (different box configurations)
12 point_num - number of evaluation points
13 max_deg_x - the maximum degree in x-direction
14 max_deg_y - the maximum degree in y-direction
15 num_knots_x - maximum number of inner knots in x direction- the real value i generated by random
16 num_knots_x - maximum number of inner knots in y direction- the real value i generated by random
17 rand_knot_num - maximum number of inserted knots- increase the multiplicity of the knots- does not add new knot values
18 max_level - maximum number of hierarchical levels
19 box_num - maximum number of boxes inserted into one level of hierarchy
20
21 This file is part of the G+Smo library.
22
23 This Source Code Form is subject to the terms of the Mozilla Public
24 License, v. 2.0. If a copy of the MPL was not distributed with this
25 file, You can obtain one at http://mozilla.org/MPL/2.0/.
26
27 Author(s):
28
29 **/
30
31 #include "gismo_unittest.h"
32
33 using namespace gismo;
34
35
36 //lvl- max level of domain refinement, max number of boxes inserted
37 //per level, basis, alignement with the prewious level
random_refinement(int lvl,int max_nb,gismo::gsHTensorBasis<2> * bas,bool aligned=false)38 std::vector<index_t> random_refinement(int lvl, int max_nb,
39 gismo::gsHTensorBasis<2> *bas,
40 bool aligned = false)
41 {
42
43 //int max_x, max_y;
44 int boxes_in_lvl, bi,span_size;
45 std::vector<index_t> q;
46 std::vector<gsVector<index_t> > boxes;// boxes stores the boxes inserted to the previous level
47 std::vector<gsVector<index_t> > boxes_new;
48 gsVector<index_t, 2> i1;
49 gsVector<index_t, 2> i2;
50 span_size = 1 << lvl;
51 //gsDebug<<"Spansize is :"<< span_size<<"\n";
52
53 // left bottom corner box
54 i1.setZero();
55 i2[0] = bas->getBases()[0]->component(0).knots().unique().size()-1;
56 i2[1] = bas->getBases()[0]->component(1).knots().unique().size()-1;
57
58 i1[0] = i1[0] << (lvl); // get indices in lvl
59 i1[1] = i1[1] << (lvl);
60 i2[0] = i2[0] << (lvl);
61 i2[1] = i2[1] << (lvl);
62
63 boxes.push_back(i1);//insert the root of the tree into the boxes
64 boxes.push_back(i2);
65 srand((unsigned)time(NULL));//seed the random alg.
66 for(int i = 0; i < lvl; i++){//insert lvl levels
67 //gsDebug<<"level: "<<i+1<< "\n";
68 boxes_in_lvl = (rand()%max_nb);
69 if(boxes_in_lvl == 0){
70 boxes_in_lvl++;
71 }
72 for(int j = 0; j < boxes_in_lvl; j++){//insert number of boxes
73 bi = rand() % (boxes.size()/2);//pick a box from prewious level
74 if(aligned){
75 //test if the box is a box not a point- insert only boxes
76 i1[0] = ( ( (rand() % (boxes[2*bi+1][0]-boxes[2*bi][0]) ) + boxes[2*bi][0]) /span_size) << i;// left bottom corner box
77 i1[1] = ( ( (rand() % (boxes[2*bi+1][1]-boxes[2*bi][1]) ) + boxes[2*bi][1]) /span_size) << i;
78 i2[0] = ( ( (rand() % (boxes[2*bi+1][0]-i1[0]) +1) + i1[0]) /span_size) << i;// right top corner of the box
79 i2[1] = ( ( (rand() % (boxes[2*bi+1][1]-i1[1]) +1) + i1[1]) /span_size) << i;
80
81 i2[0] = math::max( i2[0], i1[0] + (bas->degree(1)+2)*(1 << lvl) );
82 i2[1] = math::max( i2[1], i1[1] + (bas->degree(1)+2)*(1 << lvl) );
83
84
85 //gsDebug<<"\naligned box inserted ["<< i1[0]<<" , "<< i1[1]<<"] ["<< i2[0]<<" , "<< i2[1]<<"]"<< " to level "<<i+1 <<"\n";
86 boxes_new.push_back(i1);
87 boxes_new.push_back(i2);
88 q.push_back(i+1);
89 i1[0] = i1[0] >> (lvl - (i+1)); // get indices in level i+1
90 i1[1] = i1[1] >> (lvl - (i+1));
91 i2[0] = i2[0] >> (lvl - (i+1));
92 i2[1] = i2[1] >> (lvl - (i+1));
93 q.push_back(i1[0]);
94 q.push_back(i1[1]);
95 q.push_back(i2[0]);
96 q.push_back(i2[1]);
97 //bas->insert_box(i1,i2,i+1);
98 }else{
99 i1[0] = (rand() % (boxes[2*bi+1][0]-boxes[2*bi][0]) ) + boxes[2*bi][0];
100 i1[1] = (rand() % (boxes[2*bi+1][1]-boxes[2*bi][1]) ) + boxes[2*bi][1];
101 i2[0] = (rand() % (boxes[2*bi+1][0]-i1[0]) +1) + i1[0];
102 i2[1] = (rand() % (boxes[2*bi+1][1]-i1[1]) +1) + i1[1];
103 if( (i1[0]==i2[0]) || (i1[1]==i2[1]) ){
104 i2[0] += span_size;
105 i2[1] += span_size;
106 }
107 //gsDebug<<"\nbox inserted ["<< i1[0]<<" , "<< i1[1]<<"] ["<< i2[0]<<" , "<< i2[1]<<"]"<< " to level "<<i+1 <<"\n";
108 boxes_new.push_back(i1);
109 boxes_new.push_back(i2);
110 q.push_back(i+1);
111 i1[0] = i1[0] >> (lvl - (i+1)); // get indices in level i+1
112 i1[1] = i1[1] >> (lvl - (i+1));
113 i2[0] = i2[0] >> (lvl - (i+1));
114 i2[1] = i2[1] >> (lvl - (i+1));
115 q.push_back(i1[0]);
116 q.push_back(i1[1]);
117 q.push_back(i2[0]);
118 q.push_back(i2[1]);
119 //bas->insert_box(i1,i2,i+1);
120 }
121 }
122 span_size = span_size/2;
123 //gsDebug<<"Spansize is :"<< span_size<<endl;
124 boxes.clear();
125 for(unsigned int k = 0; k < boxes_new.size();k++){
126 boxes.push_back(boxes_new[k]);
127 }
128 boxes_new.clear();
129 }
130 return q;
131 }
132
133
SUITE(gsThbs_geometry_test)134 SUITE(gsThbs_geometry_test)
135 {
136
137 TEST(gsThbs_geometry_test)
138 {
139 gsVector<index_t> iv1;
140 gsVector<index_t> iv2;
141 iv1.resize(2);
142 iv2.resize(2);
143
144 srand((unsigned)time(NULL));//seed the random alg.
145 int test_num = 1;
146 int point_num = 10;
147 int max_deg_x = 2;
148 int max_deg_y = 2;
149 int num_knots_x = 2;
150 int num_knots_y = 2;
151 int rand_knot_num = 0;
152 int max_level = 2;
153 int box_num = 5;
154 for(int i =0; i < test_num;i++){
155 int deg_x = (rand()%max_deg_x)+1;
156 int deg_y = (rand()%max_deg_y)+1;
157 int kn = (rand()%num_knots_x)+1;
158 gsKnotVector<> T_KV (0, 1, kn , deg_x+1, 1 ) ;
159 kn = (rand()%num_knots_y)+1;
160 gsKnotVector<> T_KV1 (0, 1,kn,deg_y+1,1) ;
161 //gsInfo<<"Knot Vector"<<T_KV<<"\n";
162 //gsInfo<<"Knot Vector"<<T_KV1<<"\n";
163 for (int i2 = 0; i2 < rand_knot_num; i2 ++){
164 int kn2 = rand()%T_KV.size();
165 if(T_KV.multiplicity(T_KV[kn2]) < deg_x ){
166 T_KV.insert(T_KV[kn2]);
167 }
168 }
169 for (int i3 = 0; i3 < rand_knot_num; i3 ++){
170 int kn3 = rand()%T_KV1.size();
171 if(T_KV1.multiplicity(T_KV1[kn3]) < deg_y ){
172 T_KV1.insert(T_KV1[kn3]);
173 }
174 }
175 //gsInfo<<"Knot Vector"<<T_KV<<"\n";
176 //gsInfo<<"Knot Vector"<<T_KV1<<"\n";
177 gsTensorBSplineBasis<2> T_tbasis( T_KV, T_KV1 );
178
179 int T_nlevels =(rand()%max_level)+2;
180 gsTHBSplineBasis<2> TT ( T_tbasis ); //necessary to create the refinement
181 gsTHBSplineBasis<2> THB ( T_tbasis, random_refinement(T_nlevels-1, box_num, &TT )) ;
182
183 //gsInfo<<"knot vector"<< THB.m_cvs[0][0]<<"\n";
184 //gsInfo<<"knot vector"<< THB.m_cvs[1][0]<<"\n";
185 //gsInfo<<"Basis has degree "<< deg_x <<" and "<< deg_y <<". The number of levels is "<< T_nlevels<<"\n";
186
187 //gsInfo<<"The tree has "<< THB.tree().size() << " nodes.\n" << "\n";
188
189 // THB.printCharMatrix();
190
191 gsMatrix<> T_para = THB.support();
192 gsVector<> T_c0 = T_para.col(0);
193 gsVector<> T_c1 = T_para.col(1);
194 gsMatrix<> T_pts = uniformPointGrid(T_c0,T_c1, point_num) ;
195
196 //gsInfo<<"\n"<< "The parameter range is: "<< "\n" << T_para <<"\n";
197 ////////coefficient matrix creation///////////
198 gsMatrix<> anch = THB.anchors();
199
200 //gsInfo<<"anchors"<< anch<<"\n";
201 gsMatrix<> eye = gsMatrix<>::Identity(THB.size(), THB.size());
202 gsTHBSpline<2> THB_geometry ( THB, eye);
203 gsMatrix<> T_ev_geom = THB_geometry.eval(anch);
204 gsMatrix<> T_ev = THB.eval(anch) ;
205 //gsInfo<<"eval geometry \n"<< T_ev_geom <<"\n";
206 CHECK_MATRIX_PARTITION_OF_UNIT_CLOSE(T_ev, (real_t)0.0000001);
207 }
208
209 }
210
211 }
212