G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsTHBSpline.hpp
Go to the documentation of this file.
1
15#pragma once
16
17#include <gsNurbs/gsBoehm.h>
18
20
21namespace gismo
22{
23
24// ************************************************
25// Public member functions
26// ************************************************
27
28template<short_t d, class T>
30{
31 typedef typename gsHDomain<d>::point point;
32
33 const gsHDomain<d>& tree = this->basis().tree();
34
35 // Construct a box covering the whole parameter domain.
36 const point & uCornerGlob = tree.upperCorner();
37 point uCornerLoc;
38
39 index_t maxInsLevel = tree.getMaxInsLevel();
40 tree.global2localIndex(uCornerGlob, maxInsLevel, uCornerLoc);
41
42 std::vector<index_t> wholeDomainAsBox(2*d+1, 0);
43 wholeDomainAsBox[0] = maxInsLevel;
44
45 std::copy(uCornerLoc.data(), uCornerLoc.data()+d, wholeDomainAsBox.begin()+d+1);
46
47 // Refine the whole domain to the finest level present there.
48 this->refineElements( wholeDomainAsBox );
49
50 tensorBasis & tpBasis = this->basis().tensorLevel(this->basis().maxLevel());
51
52 // makeGeometry returns an abstract class, therefore we need to cast to the particular.
53 result = give( static_cast<gsTensorBSpline<d,T>&>(*tpBasis.makeGeometry(this->coefs())) );
54}
55
56template<short_t d, class T>
57void gsTHBSpline<d, T>::increaseMultiplicity(index_t lvl, int dir, T knotValue, int mult)
58{
59 gsWarn<<"gsTHBSpline<d, T>::increaseMultiplicity: This code is not working properly!"<<std::endl;
60 // Copy the current characteristic matrices
61 gsDebug<<"in geo"<<std::endl;
62 std::vector<gsSortedVector<index_t> > OX = this->basis().getXmatrix();
63
64 // Insert the knot in the basis
65 this->basis().increaseMultiplicity(lvl,dir,knotValue,mult);
66 gsDebug<<"increased"<<std::endl;
67 // Compute the transfer matrix
68 gsSparseMatrix<T> trMatrix;
69 this->basis().transfer(OX, trMatrix);
70 gsDebug<<"transfer"<<std::endl;
71 // Multiply the coeffs by the transfer matrix
72 this->m_coefs = trMatrix * this->m_coefs;
73}
74
75// Return the list of B-spline patches to represent a THB-spline geometry
76// \returns b1 bottom left corners of the box (vector of indices with respect to the gsCompactKnotVector of the highest possible level)
77// \returns b2 top right corners of the box (vector of indices with respect to the gsCompactKnotVector of the highest possible level)
78// \returns level levels of the boxes (level[i]: level of the i-th box,)
79// \returns bpatches list of B-spline patches associated with the boxes
80template<short_t d, class T>
81//void gsTHBSplineBasis<d,T>::getBsplinePatches(gsMatrix<T>& geom_coef, gsMatrix<T>& cp, gsMatrix<index_t>& b1, gsMatrix<index_t>& b2, gsVector<index_t>& level, gsMatrix<index_t>& nvertices) const
82//void gsTHBSpline<d,T>::getBsplinePatches(gsMatrix<index_t>& b1, gsMatrix<index_t>& b2, gsVector<index_t>& level, std::vector< gsTensorBSpline<2> >& bpatches) const
84{
85 GISMO_UNUSED(b1);
86 GISMO_UNUSED(b2);
87 GISMO_UNUSED(level);
89 //------------------------------------------------------------------------------------------------------------------------------
90 // define the 6 boxes and the corresponding levels for the MTU fillet example *** TO BE SUBSTITUTED WITH AUTOMATIC SPLITTING ***
91 //------------------------------------------------------------------------------------------------------------------------------
92 /*int nboxes = 6;
93 level.resize(nboxes);
94 b1.resize(nboxes,this->dim());
95 b2.resize(nboxes,this->dim());
96
97 level[0] = 0; b1(0,0) = 0; b1(0,1) = 0; b2(0,0) = 16; b2(0,1) = 80;
98 level[1] = 1; b1(1,0) = 16; b1(1,1) = 0; b2(1,0) = 40; b2(1,1) = 80;
99 level[2] = 2; b1(2,0) = 40; b1(2,1) = 0; b2(2,0) = 60; b2(2,1) = 80;
100 level[3] = 2; b1(3,0) = 60; b1(3,1) = 0; b2(3,0) = 80; b2(3,1) = 16;
101 level[4] = 2; b1(4,0) = 60; b1(4,1) = 66; b2(4,0) = 80; b2(4,1) = 80;
102 level[5] = 3; b1(5,0) = 60; b1(5,1) = 16; b2(5,0) = 80; b2(5,1) = 66;*/
103 //------------------------------------------------------------------------------------------------------------------------------
104 // define the 8 boxes and the corresponding levels for the face example *** TO BE SUBSTITUTED WITH AUTOMATIC SPLITTING ***
105 //------------------------------------------------------------------------------------------------------------------------------
106 /*int nboxes = 8;
107 level.resize(nboxes);
108 b1.resize(nboxes,this->dim());
109 b2.resize(nboxes,this->dim());
110
111 level[0] = 1; b1(0,0) = 0; b1(0,1) = 20; b2(0,0) = 32; b2(0,1) = 24;
112 level[1] = 2; b1(1,0) = 0; b1(1,1) = 0; b2(1,0) = 48; b2(1,1) = 20;
113 level[2] = 2; b1(2,0) = 32; b1(2,1) = 20; b2(2,0) = 48; b2(2,1) = 24;
114 level[3] = 2; b1(3,0) = 0; b1(3,1) = 24; b2(3,0) = 12; b2(3,1) = 96;
115 level[4] = 2; b1(4,0) = 12; b1(4,1) = 24; b2(4,0) = 36; b2(4,1) = 48;
116 level[5] = 2; b1(5,0) = 36; b1(5,1) = 24; b2(5,0) = 48; b2(5,1) = 96;
117 level[6] = 2; b1(6,0) = 12; b1(6,1) = 82; b2(6,0) = 36; b2(6,1) = 96;
118 level[7] = 3; b1(7,0) = 12; b1(7,1) = 48; b2(7,0) = 36; b2(7,1) = 82;*/
119 //------------------------------------------------------------------------------------------------------------------------------
120 // define the 5 boxes and the corresponding levels for file 15 *** TO BE SUBSTITUTED WITH AUTOMATIC SPLITTING ***
121 //------------------------------------------------------------------------------------------------------------------------------
122 /*int nboxes = 5;
123 level.resize(nboxes);
124 b1.resize(nboxes,this->dim());
125 b2.resize(nboxes,this->dim());
126
127 level[0] = 0; b1(0,0) = 0; b1(0,1) = 0; b2(0,0) = 16; b2(0,1) = 24; //int l = 0; b1[0] = 0; b1[1] = 0; b2[0] = 16; b2[1] = 24;
128
129 level[1] = 1; b1(1,0) = 16; b1(1,1) = 0; b2(1,0) = 28; b2(1,1) = 8; //int l = 1; b1[0] = 16; b1[1] = 0; b2[0] = 28; b2[1] = 8;
130 level[2] = 1; b1(2,0) = 16; b1(2,1) = 16; b2(2,0) = 28; b2(2,1) = 24; //int l = 1; b1[0] = 16; b1[1] = 16; b2[0] = 28; b2[1] = 24;
131 level[3] = 1; b1(3,0) = 24; b1(3,1) = 8; b2(3,0) = 28; b2(3,1) = 16; //int l = 1; b1[0] = 24; b1[1] = 0; b2[0] = 28; b2[1] = 16;
132
133 level[4] = 2; b1(4,0) = 16; b1(4,1) = 8; b2(4,0) = 24; b2(4,1) = 16;*/ //int l = 2; b1[0] = 16; b1[1] = 8; b2[0] = 24; b2[1] = 16;
134 //------------------------------------------------------------------------------------------------------------------------------
135 // AUTOMATIC SPLITTING
136 //------------------------------------------------------------------------------------------------------------------------------
137 //this->basis().component(1); //m_tree.getBoxes(b1,b2,level);
138 //this->basis().component(i).degree();
139 //int nboxes = level.size();
140 //------------------------------------------------------------------------------------------------------------------------------
141 // iteration on the boxes to call getBsplinePatchGlobal()
142 //------------------------------------------------------------------------------------------------------------------------------
143 /*gsVector<index_t> p1, p2;
144 p1.resize(this->dim());
145 p2.resize(this->dim());
146 gsMatrix<T> temp1, temp2;
147 gsCompactKnotVector<T> cku, ckv;
148 nvertices.resize(nboxes,this->dim());
149
150 for (int i = 0; i < nboxes; i++){
151 p1(0) = b1(i,0); p1(1) = b1(i,1); p2(0) = b2(i,0); p2(1) = b2(i,1);
152 this->getBsplinePatchGlobal(p1, p2, level[i], geom_coef, temp1, cku, ckv);
153
154 gsDebug<<"cku: "<<cku<<std::endl<<"ckv: "<<ckv<<std::endl;
155
156 if (i == 0){
157 cp = temp1;
158 }else{
159 int cprows = cp.rows();
160 temp2.resize(cp.rows()+temp1.rows(),cp.cols()); // is there a way to append the temp1 matrix to cp without temp2?
161
162 for (int j = 0; j < cp.rows(); j++){
163 for (int k = 0; k < cp.cols(); k++){
164 temp2(j,k) = cp(j,k);
165 }
166 }
167
168 for (int j = 0; j < temp1.rows(); j++){
169 for (int k = 0; k < temp1.cols(); k++){
170 //gsDebug<<cprows+j<<endl;
171 temp2(cprows+j,k) = temp1(j,k);
172 }
173 }
174 cp = temp2;
175 }
176
177 nvertices(i,0) = cku.size()-cku.degree()-1;
178 nvertices(i,1) = ckv.size()-ckv.degree()-1;
179
180 }*/
181 //gsDebug<<cp<<endl<<endl;
182 //gsDebug<<cp.rows()<<endl;//ok
183 //gsDebug<<"b1 "<<b1.rows()<<endl;
184 //gsDebug<<"b2 "<<b2.rows()<<endl;
185 //gsDebug<<"level "<<level.rows()<<endl;
186 //gsDebug<<"nv "<<nvertices.rows()<<endl;//ko
187}
188// ************************************************
189// Private member functions
190// ************************************************
191
192//
193// Return the B-spline representation of a THB-spline subpatch
194// \param b1 bottom left corner of the box (vector of indices with respect to the gsCompactKnotVector of the highest possible level)
195// \param b2 top right corner of the box (vector of indices with respect to the gsCompactKnotVector of the highest possible level)
196// \param level level of the box
197// \param geom_coef control points of the THB-spline geometry
198// \returns cp control control points of the B-spline patch
199// \returns k1 knot vector of the B-spline patch (first dimension)
200// \returns k2 knot vector of the B-spline patch (second dimension)
201/*template<short_t d, class T>
202void gsTHBSplineBasis<d,T>::getBsplinePatchGlobal(gsVector<index_t> b1, gsVector<index_t> b2, unsigned level, gsMatrix<T>& geom_coef, gsMatrix<T>& cp, gsCompactKnotVector<T>& k1, gsCompactKnotVector<T>& k2) const
203{
204 // check if the indices in b1, and b2 are correct with respect to the given level
205 if(b1[0]%Qlocal2global(1,0, level) != 0 ){
206 b1[0] -= b1[0]%Qlocal2global(1,0, level);
207 }
208 if(b1[1]%Qlocal2global(1,0, level) != 0 ){
209 b1[1] -= (b1[1]%Qlocal2global(1,0, level));
210 }
211 if(b2[0]%Qlocal2global(1,0, level) != 0){
212 b2[0] += Qlocal2global(1,0, level) -(b2[0]%Qlocal2global(1,0, level));
213 }
214 if(b2[1]%Qlocal2global(1,0, level) != 0){
215 b2[1] += Qlocal2global(1,0, level) -(b2[1]%Qlocal2global(1,0, level));
216 }
217 // select the indices of all B-splines of the given level acting on the given box
218 int i0 = Qglobal2local(b1[0],level,this->m_tree.m_index_level);
219 i0 = this->m_bases[level]->component(0).knots().lastKnotIndex(i0) - this->m_deg[0];
220 int i1 = Qglobal2local(b2[0],level,this->m_tree.m_index_level);
221 i1 = this->m_bases[level]->component(0).knots().firstKnotIndex(i1) - 1;
222
223 int j0 = Qglobal2local(b1[1],level,this->m_tree.m_index_level);
224 j0 = this->m_bases[level]->component(1).knots().lastKnotIndex(j0) - this->m_deg[1];
225 int j1 = Qglobal2local(b2[1],level,this->m_tree.m_index_level);
226 j1 = this->m_bases[level]->component(1).knots().firstKnotIndex(j1) - 1;
227
228 for(int i = 0; i < geom_coef.cols(); i++){
229 initialize_cmatrix(geom_coef, i, level);
230 gsMatrix<T> *temp = new gsMatrix<T>(1,1);
231 globalRefinement(level, *temp);
232
233 for(int k = i0; k <= i1; k++){
234 for(int j = j0; j <= j1; j++){
235 (*temp)(j-j0,k-i0) = (*temp)(j,k);
236 }
237 }
238 temp->conservativeResize(j1-j0+1,i1-i0+1);
239
240 if(i == 0){
241 cp.resize(temp->cols()*temp->rows(), geom_coef.cols());
242 }
243 return_cp_1D(*temp, i, cp);
244 }
245 // compute the new vectors for the B-spline patch
246 k1 = gsCompactKnotVector<T>(this->m_deg[0], this->m_bases[level]->component(0).knots().begin() + i0 , this->m_bases[level]->component(0).knots().begin() + i1 + this->m_deg[0] + 2);
247 k2 = gsCompactKnotVector<T>(this->m_deg[1], this->m_bases[level]->component(1).knots().begin() + j0 , this->m_bases[level]->component(1).knots().begin() + j1 + this->m_deg[1] + 2);
248}*/
249
250//
251// Function called by getBsplinePatchGlobal
252// \param level
253/*template<short_t d, class T>
254void gsTHBSplineBasis<d,T>::globalRefinement(int level, gsMatrix<T>& coeffs)const
255{
256 //coeffs.resize(this->m_bases[0]->component(1).knots().size()-this->m_deg[1]-1,this->m_bases[0]->component(0).knots().size()-this->m_deg[0]-1);
257 coeffs.resize(this->m_bases[0]->size(1),this->m_bases[0]->size(0));
258 for(int j = 0; j < coeffs.rows(); j++){
259 for(int k = 0; k < coeffs.cols(); k++){
260 unsigned s = this->m_bases[0]->index(k,j);//this->fromPair(k,j, 0);
261 if(this->m_xmatrix[0].count( s ) > 0){
262 coeffs(j,k) = this->m_cmatrix[0].find( s )->second;
263 }else{
264 coeffs(j,k) = 0;
265 }
266 }
267 }
268
269 for(int l = 1; l <=level; l++){
270
271 // global dyadic refinement with respect to previous level
272 gsCompactKnotVector<T>k1;
273 gsCompactKnotVector<T>k2;
274 k1 = this->m_bases[l-1]->component(0).knots();
275 k2 = this->m_bases[l-1]->component(1).knots();
276 std::vector<T> knots_x;
277 std::vector<T> knots_y;
278 for(unsigned int i = 1; i < this->m_bases[l]->component(0).knots().unique().size(); i = i+2){
279 knots_x.push_back(this->m_bases[l]->component(0).knots().unique()[i]);
280 }
281 for(unsigned int i = 1; i < this->m_bases[l]->component(1).knots().unique().size(); i = i+2){
282 knots_y.push_back(this->m_bases[l]->component(1).knots().unique()[i]);
283 }
284
285 gsBoehmRefine(k2,coeffs,this->m_deg[1],knots_y);
286 coeffs.transposeInPlace();
287
288 gsBoehmRefine(k1,coeffs,this->m_deg[0],knots_x);
289 coeffs.transposeInPlace();
290
291 //overwrite the whole matrix
292 for(int j = 0; j < coeffs.rows(); j++){
293 for(int k = 0; k < coeffs.cols(); k++){
294 unsigned s = this->m_bases[l]->index(k,j);
295 if(this->m_xmatrix[l].count( s ) > 0){
296 coeffs(j,k) = this->m_cmatrix[l].find( s )->second;
297 }
298 }
299 }
300 }
301}*/
302
303//
304// Initializes the m_cmatrix up to given level with the coeffs of the geometry
305// \param col dimension (0, 1, 2 = x, y, z)
306// \param c_level level
307/*template<short_t d, class T>
308void gsTHBSplineBasis<d,T>::initialize_cmatrix(gsMatrix<T>&geom_coeff, int col, int c_level) const{
309 int counter = 0;
310 for(int i = 0; i <= c_level; i++){
311 this->m_cmatrix.push_back( std::map <index_t, T>());
312 for(typename CMatrix::const_iterator it =
313 this->m_xmatrix[i].begin(); it != this->m_xmatrix[i].end(); it++){
314 this->m_cmatrix[i][ *it ] = geom_coeff(counter,col);
315 counter++;
316 }
317 }
318}*/
319
320//
321// Converts the coefficient matrix to a column of the control points matrix with respect to the given direction
322// \param mat
323// \param direction dimension (0, 1, 2 = x, y, z)
324// \returns cp the column direction of cp is updated with the related coordinate of the control points
325/*template<short_t d, class T>
326void gsTHBSplineBasis<d,T>::return_cp_1D(const gsMatrix<T> & mat, int direction, gsMatrix<T>& cp)const{
327 GISMO_ASSERT((mat.cols()*mat.rows() == cp.rows()), "Wrong matrix dimension.");
328 int counter = 0;
329 for(int j = 0; j < mat.rows(); j++){
330 for(int i = 0; i < mat.cols(); i++){
331 cp(counter, direction) = mat(j,i);
332 counter ++;
333 }
334 }
335}*/
336
337template<short_t d, class T>
338void gsTHBSpline<d,T>::slice(index_t dir_fixed,T par,
339 typename gsTHBSpline<d,T>::BoundaryGeometryType & result) const
340{
341 GISMO_ASSERT(d-1>=0,"d must be greater or equal than 1");
342 GISMO_ASSERT(dir_fixed>=0 && static_cast<index_t>(dir_fixed)<d,"cannot fix a dir greater than dim or smaller than 0");
343
344 typedef typename gsTHBSpline<d,T>::BoundaryBasisType THBBoundaryBasis;
345 typedef typename gsTHBSpline<d,T>::BoundaryGeometryType THBBoundary;
346
347 const THBBoundaryBasis * bBasis = this->basis().basisSlice(dir_fixed,par);
348 if(d==1)
349 {
350 gsMatrix<T> val(1,1),point;
351 val(0,0)=par;
352 this->eval_into(val,point);
353 result = THBBoundary(*bBasis,point);
354 }
355 else
356 {
357 gsMatrix<T> vals,anchorsSlice,anchorsInGeom;
358 bBasis->anchors_into(anchorsSlice);
359 anchorsInGeom.resize(anchorsSlice.rows()+1,anchorsSlice.cols());
360 anchorsInGeom.topRows(dir_fixed)=anchorsSlice.topRows(dir_fixed);
361 anchorsInGeom.row(dir_fixed)=gsVector<T>::Constant(anchorsSlice.cols(),par);
362 anchorsInGeom.bottomRows(anchorsSlice.rows()-dir_fixed)=anchorsSlice.bottomRows(anchorsSlice.rows()-dir_fixed);
363 this->eval_into(anchorsInGeom,vals);
364 THBBoundary* geom =
365 dynamic_cast<THBBoundary*>(bBasis->interpolateData(vals,anchorsSlice).release());
366 GISMO_ASSERT(geom!=NULL,"bBasis should have BoundaryGeometryType.");
367 result = *geom;
368 delete geom;
369 }
370 delete bBasis;
371}
372
373
374namespace internal
375{
377template<short_t d, class T>
378class gsXml< gsTHBSpline<d,T> >
379{
380private:
381 gsXml() { }
382public:
383 GSXML_COMMON_FUNCTIONS(gsTHBSpline<TMPLA2(d,T)>);
384 static std::string tag () { return "Geometry"; }
385 static std::string type () { return "THBSpline"+to_string(d); }
386
387 static gsTHBSpline<d,T> * get (gsXmlNode * node)
388 {
389 return getGeometryFromXml< gsTHBSpline<d,T> >(node);
390 }
391
392 static gsXmlNode * put (const gsTHBSpline<d,T> & obj,
393 gsXmlTree & data )
394 {
395 return putGeometryToXml< gsTHBSpline<d,T> >(obj,data);
396 }
397};
398
399}// namespace internal
400
401}// namespace gismo
virtual memory::unique_ptr< gsGeometry< T > > makeGeometry(gsMatrix< T > coefs) const =0
Create a gsGeometry of proper type for this basis with the given coefficient matrix.
Class with a hierarchical domain structure represented by a box k-d-tree.
Definition gsHDomain.h:76
const point & upperCorner() const
Accessor for gsHDomain::m_upperIndex.
Definition gsHDomain.h:249
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
Sparse matrix class, based on gsEigen::SparseMatrix.
Definition gsSparseMatrix.h:139
A truncated hierarchical B-Spline function, in d dimensions.
Definition gsTHBSpline.h:38
void convertToBSpline(gsTensorBSpline< d, T > &result)
Refines the whole domain to the finest level present in the mesh. Returns the refined geometry as res...
Definition gsTHBSpline.hpp:29
void getBsplinePatches(gsMatrix< index_t > &b1, gsMatrix< index_t > &b2, gsVector< index_t > &level) const
Return the list of B-spline patches to represent a THB-spline geometry.
Definition gsTHBSpline.hpp:83
void increaseMultiplicity(index_t lvl, int dir, T knotValue, int mult=1)
Increase multiplicity of knot-value knotValue in level lvl and direction dir by mult.
Definition gsTHBSpline.hpp:57
void slice(index_t dir_fixed, T par, BoundaryGeometryType &result) const
Definition gsTHBSpline.hpp:338
A tensor product B-spline basis.
Definition gsTensorBSplineBasis.h:37
A tensor product of d B-spline functions, with arbitrary target dimension.
Definition gsTensorBSpline.h:45
Boehm's algorithm for knot insertion.
#define index_t
Definition gsConfig.h:32
Provides declaration of ConstantFunction class.
#define gsDebug
Definition gsDebug.h:61
#define GISMO_NO_IMPLEMENTATION
Definition gsDebug.h:129
#define gsWarn
Definition gsDebug.h:50
#define GISMO_UNUSED(x)
Definition gsDebug.h:112
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
std::string to_string(const unsigned &i)
Helper to convert small unsigned to string.
Definition gsXml.cpp:74
The G+Smo namespace, containing all definitions for the library.
S give(S &x)
Definition gsMemory.h:266