NURBS3DSurface.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2007-2019 PCOpt/NTUA
9  Copyright (C) 2013-2019 FOSS GP
10  Copyright (C) 2019-2020 OpenCFD Ltd.
11 -------------------------------------------------------------------------------
12 License
13  This file is part of OpenFOAM.
14 
15  OpenFOAM is free software: you can redistribute it and/or modify it
16  under the terms of the GNU General Public License as published by
17  the Free Software Foundation, either version 3 of the License, or
18  (at your option) any later version.
19 
20  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23  for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include "NURBS3DSurface.H"
31 #include "vtkSurfaceWriter.H"
32 #include "OFstream.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 label NURBS3DSurface::sgn(const scalar val) const
42 {
43  return val >= scalar(0) ? 1 : -1;
44 }
45 
46 
47 scalar NURBS3DSurface::abs(const scalar val) const
48 {
49  return (sgn(val) == 1)? val: -val;
50 }
51 
52 
53 label NURBS3DSurface::mod(const label x, const label interval) const
54 {
55  label ratio(x%interval);
56  return ratio < 0 ? ratio+interval : ratio;
57 }
58 
59 
60 void NURBS3DSurface::setCPUVLinking()
61 {
62  const label uNCPs(uBasis_.nCPs());
63  const label vNCPs(vBasis_.nCPs());
64 
65  CPsUCPIs_.setSize(uNCPs*vNCPs, -1);
66  CPsVCPIs_.setSize(uNCPs*vNCPs, -1);
67 
68  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
69  {
70  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
71  {
72  const label CPI(vCPI*uNCPs + uCPI);
73  CPsUCPIs_[CPI] = uCPI;
74  CPsVCPIs_[CPI] = vCPI;
75  }
76  }
77 }
78 
79 
80 void NURBS3DSurface::setUniformUV
81 (
82  scalarList& u,
83  scalarList& v,
84  const label nUPts,
85  const label nVPts
86 ) const
87 {
88  u.setSize(nUPts*nVPts, Zero);
89  v.setSize(nUPts*nVPts, Zero);
90 
91  for (label uI = 0; uI<nUPts; uI++)
92  {
93  scalar uValue = scalar(uI)/scalar(nUPts - 1);
94  for (label vI = 0; vI<nVPts; vI++)
95  {
96  const label ptI(uI*nVPts + vI);
97  u[ptI] = uValue;
98  v[ptI] = scalar(vI)/scalar(nVPts - 1);
99  }
100  }
101 }
102 
103 
104 void NURBS3DSurface::setUniformUV()
105 {
106  setUniformUV(u_, v_, nUPts_, nVPts_);
107 }
108 
109 
110 bool NURBS3DSurface::bound
111 (
112  scalar& u,
113  scalar& v,
114  const scalar minVal,
115  const scalar maxVal
116 ) const
117 {
118  bool boundPoint =
119  boundDirection(u, minVal, maxVal)
120  || boundDirection(v, minVal, maxVal);
121 
122  return boundPoint;
123 }
124 
125 
126 bool NURBS3DSurface::boundDirection
127 (
128  scalar& u,
129  const scalar minVal,
130  const scalar maxVal
131 ) const
132 {
133  bool boundPoint(false);
134 
135  if (u < scalar(0))
136  {
137  u = minVal;
138  boundPoint = true;
139  }
140  else if (u > scalar(1))
141  {
142  u = maxVal;
143  boundPoint = true;
144  }
145 
146  return boundPoint;
147 }
148 
149 
150 void NURBS3DSurface::setEquidistantR
151 (
152  scalarList& R,
153  const scalar SHeld,
154  const label paramR,
155  const label lenAcc = 25,
156  const label maxIter = 10,
157  const label spacingCorrInterval = -1,
158  const scalar tolerance = 1.e-5
159 ) const
160 {
161  const label nPts(R.size());
162  scalar SNull(SHeld);
163  scalar xLength(Zero);
164  const scalar rLength(scalar(1) / scalar(nPts - 1));
165 
166  if (paramR == PARAMU)
167  {
168  xLength = lengthU(SHeld) / (nPts - 1);
169  }
170  else
171  {
172  xLength = lengthV(SHeld) / (nPts - 1);
173  }
174 
175  R[0] = Zero;
176  R[nPts-1] = scalar(1);
177 
178  for (label ptI=1; ptI<(nPts - 1); ptI++)
179  {
180  const scalar& RPrev(R[ptI - 1]);
181  scalar& RCurr(R[ptI]);
182  scalar direc(1);
183  scalar xDiff(0);
184  scalar delta(0);
185  bool overReached(false);
186 
187  RCurr = RPrev + rLength;
188 
189  // Find the starting U value to ensure target is within 1 rLength.
190  while (true)
191  {
192  bool bounded(false);
193 
194  if (paramR == PARAMU)
195  {
196  bounded = bound(RCurr, SNull);
197  delta = lengthU(SHeld, RPrev, RCurr, lenAcc);
198  }
199  else
200  {
201  bounded = bound(SNull, RCurr);
202  delta = lengthV(SHeld, RPrev, RCurr, lenAcc);
203  }
204 
205  xDiff = xLength - delta;
206 
207  // Found the point.
208  if (abs(xDiff) < tolerance)
209  {
210  break;
211  }
212  else
213  {
214  direc = sgn(xDiff);
215 
216  if (bounded && (direc == 1))
217  {
218  // rLength addition makes U exceed 1 so it becomes bounded.
219  // However, the desired x location still exceeds how far the
220  // bounded rLength can move (~e-5 error).
221  // Must force U to be u=0.999999.
222  overReached = true;
223  break;
224  }
225  else if (direc == scalar(1))
226  {
227  RCurr += rLength;
228  }
229  else
230  {
231  break;
232  }
233  }
234  }
235 
236  if (!overReached)
237  {
238  label iter(0);
239 
240  while (iter < maxIter)
241  {
242  // Set the new search length to ensure convergence and next v.
243  direc /= scalar(2);
244  RCurr += direc * rLength;
245 
246  if (paramR == PARAMU)
247  {
248  bound(RCurr, SNull);
249  }
250  else
251  {
252  bound(SNull, RCurr);
253  }
254 
255  // Can employ an occasional tolerance check from beg of curve.
256  if
257  (
258  (spacingCorrInterval != -1)
259  && (mod(ptI, spacingCorrInterval) == 0)
260  )
261  {
262  if (paramR == PARAMU)
263  {
264  delta = lengthU(SHeld, Zero, RCurr, ptI*lenAcc);
265  }
266  else
267  {
268  delta = lengthV(SHeld, Zero, RCurr, ptI*lenAcc);
269  }
270 
271  xDiff = (ptI * xLength) - delta;
272  }
273  else
274  {
275  if (paramR == PARAMU)
276  {
277  delta = lengthU(SHeld, RPrev, RCurr, lenAcc);
278  }
279  else
280  {
281  delta = lengthV(SHeld, RPrev, RCurr, lenAcc);
282  }
283 
284  xDiff = xLength - delta;
285  }
286 
287  // Break if found point or find the new search direction.
288  if (abs(xDiff) < tolerance)
289  {
290  break;
291  }
292  else
293  {
294  const scalar oldDirec(direc);
295  direc = sgn(xDiff) * abs(oldDirec);
296  }
297 
298  ++iter;
299  }
300  }
301  }
302 }
303 
304 
305 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
306 
307 NURBS3DSurface::NURBS3DSurface
308 (
309  const List<vector>& CPs,
310  const label nUPts,
311  const label nVPts,
312  const NURBSbasis& uBasis,
313  const NURBSbasis& vBasis,
314  const word name
315 )
316 :
317  vectorField (nUPts*nVPts, Zero),
318 
319  CPs_(CPs),
320  u_(nUPts*nVPts, Zero),
321  v_(nUPts*nVPts, Zero),
322  weights_(CPs.size(), scalar(1)),
323  nUPts_(nUPts),
324  nVPts_(nVPts),
325  name_(name),
326  uBasis_(uBasis),
327  vBasis_(vBasis),
328 
329  givenInitNrm_(Zero),
330 
331  CPsUCPIs_(0),
332  CPsVCPIs_(0),
333 
334  nrmOrientation_(ALIGNED),
335 
336  boundaryCPIDs_(nullptr)
337 {
338  setUniformUV();
339  buildSurface();
340  setCPUVLinking();
341 }
342 
343 
344 NURBS3DSurface::NURBS3DSurface
345 (
346  const List<vector>& CPs,
347  const List<scalar>& weights,
348  const label nUPts,
349  const label nVPts,
350  const NURBSbasis& uBasis,
351  const NURBSbasis& vBasis,
352  const word name
353 )
354 :
355  vectorField(nUPts*nVPts, Zero),
356 
357  CPs_(CPs),
358  u_(nUPts*nVPts, Zero),
359  v_(nUPts*nVPts, Zero),
360  weights_(weights),
361  nUPts_(nUPts),
362  nVPts_(nVPts),
363  name_ (name),
364  uBasis_(uBasis),
365  vBasis_(vBasis),
366 
367  givenInitNrm_(Zero),
368 
369  CPsUCPIs_(0),
370  CPsVCPIs_(0),
371 
372  nrmOrientation_(ALIGNED),
373 
374  boundaryCPIDs_(nullptr)
375 {
376  setUniformUV();
377  buildSurface();
378  setCPUVLinking();
379 }
380 
381 
382 NURBS3DSurface::NURBS3DSurface
383 (
384  const List<vector>& CPs,
385  const label nUPts,
386  const label nVPts,
387  const label uDegree,
388  const label vDegree,
389  const label nCPsU,
390  const label nCPsV,
391  const word name
392 )
393 :
394  vectorField(nUPts*nVPts, Zero),
395 
396  CPs_(CPs),
397  u_(nUPts*nVPts, Zero),
398  v_(nUPts*nVPts, Zero),
399  weights_(CPs.size(), scalar(1)),
400  nUPts_(nUPts),
401  nVPts_(nVPts),
402  name_(name),
403  uBasis_(nCPsU, uDegree),
404  vBasis_(nCPsV, vDegree),
405 
406  givenInitNrm_(Zero),
407 
408  CPsUCPIs_(0),
409  CPsVCPIs_(0),
410 
411  nrmOrientation_(ALIGNED),
412 
413  boundaryCPIDs_(nullptr)
414 {
415  // Sanity checks
416  if (nCPsU*nCPsV != CPs_.size())
417  {
418  FatalErrorInFunction
419  << "nCPsU*nCPsV " << nCPsU*nCPsV
420  << " not equal to size of CPs " << CPs_.size()
421  << exit(FatalError);
422  }
423  // Initialize surface
424  setUniformUV();
425  buildSurface();
426  setCPUVLinking();
427 }
428 
429 
430 NURBS3DSurface::NURBS3DSurface
431 (
432  const List<vector>& CPs,
433  const label nUPts,
434  const label nVPts,
435  const label uDegree,
436  const label vDegree,
437  const label nCPsU,
438  const label nCPsV,
439  const scalarField& knotsU,
440  const scalarField& knotsV,
441  const word name
442 )
443 :
444  vectorField(nUPts*nVPts, Zero),
445 
446  CPs_(CPs),
447  u_(nUPts*nVPts, Zero),
448  v_(nUPts*nVPts, Zero),
449  weights_(CPs.size(), scalar(1)),
450  nUPts_(nUPts),
451  nVPts_(nVPts),
452  name_(name),
453  uBasis_(nCPsU, uDegree, knotsU),
454  vBasis_(nCPsV, vDegree, knotsV),
455 
456  givenInitNrm_(Zero),
457 
458  CPsUCPIs_(0),
459  CPsVCPIs_(0),
460 
461  nrmOrientation_(ALIGNED),
462 
463  boundaryCPIDs_(nullptr)
464 {
465  // Sanity checks
466  if (nCPsU*nCPsV != CPs_.size())
467  {
468  FatalErrorInFunction
469  << "nCPsU*nCPsV " << nCPsU*nCPsV
470  << " not equal to size of CPs " << CPs_.size()
471  << exit(FatalError);
472  }
473  // initialize surface
474  setUniformUV();
475  buildSurface();
476  setCPUVLinking();
477 }
478 
479 
480 NURBS3DSurface::NURBS3DSurface
481 (
482  const List<vector>& CPs,
483  const List<scalar>& weights,
484  const label nUPts,
485  const label nVPts,
486  const label uDegree,
487  const label vDegree,
488  const label nCPsU,
489  const label nCPsV,
490  const word name
491 )
492 :
493  vectorField(nUPts*nVPts, Zero),
494 
495  CPs_(CPs),
496  u_(nUPts*nVPts, Zero),
497  v_(nUPts*nVPts, Zero),
498  weights_(weights),
499  nUPts_(nUPts),
500  nVPts_(nVPts),
501  name_(name),
502  uBasis_(nCPsU, uDegree),
503  vBasis_(nCPsV, vDegree),
504 
505  givenInitNrm_(Zero),
506 
507  CPsUCPIs_(0),
508  CPsVCPIs_(0),
509 
510  nrmOrientation_(ALIGNED),
511 
512  boundaryCPIDs_(nullptr)
513 {
514  // Sanity checks
515  if (nCPsU*nCPsV != CPs_.size())
516  {
517  FatalErrorInFunction
518  << "nCPsU*nCPsV " << nCPsU*nCPsV
519  << " not equal to size of CPs " << CPs_.size()
520  << exit(FatalError);
521  }
522 
523  // Initialize surface
524  setUniformUV();
525  buildSurface();
526  setCPUVLinking();
527 }
528 
529 
530 NURBS3DSurface::NURBS3DSurface
531 (
532  const List<vector>& CPs,
533  const List<scalar>& weights,
534  const label nUPts,
535  const label nVPts,
536  const label uDegree,
537  const label vDegree,
538  const label nCPsU,
539  const label nCPsV,
540  const scalarField& knotsU,
541  const scalarField& knotsV,
542  const word name
543 )
544 :
545  vectorField(nUPts*nVPts, Zero),
546 
547  CPs_(CPs),
548  u_(nUPts*nVPts, Zero),
549  v_(nUPts*nVPts, Zero),
550  weights_(weights),
551  nUPts_(nUPts),
552  nVPts_(nVPts),
553  name_(name),
554  uBasis_(nCPsU, uDegree, knotsU),
555  vBasis_(nCPsV, vDegree, knotsV),
556 
557  givenInitNrm_(Zero),
558 
559  CPsUCPIs_(0),
560  CPsVCPIs_(0),
561 
562  nrmOrientation_(ALIGNED),
563 
564  boundaryCPIDs_(nullptr)
565 {
566  // Sanity checks
567  if (nCPsU*nCPsV != CPs_.size())
568  {
569  FatalErrorInFunction
570  << "nCPsU*nCPsV " << nCPsU*nCPsV
571  << " not equal to size of CPs " << CPs_.size()
572  << exit(FatalError);
573  }
574  // Initialize surface
575  setUniformUV();
576  buildSurface();
577  setCPUVLinking();
578 }
579 
580 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
581 
582 // * * * * * * * * * * * * * * * * Set Functions * * * * * * * * * * * * * * //
583 
584 void NURBS3DSurface::setNrmOrientation
585 (
586  const vector& givenNrm,
587  const scalar u,
588  const scalar v
589 )
590 {
591  vector surfaceNrm(surfaceDerivativeU(u, v) ^ surfaceDerivativeV(u, v));
592 
593  givenInitNrm_ = givenNrm;
594  surfaceNrm /= mag(surfaceNrm);
595 
596  const scalar relation(givenNrm & surfaceNrm);
597 
598  if (relation >= 0)
599  {
600  nrmOrientation_ = ALIGNED;
601  }
602  else
603  {
604  nrmOrientation_ = OPPOSED;
605  }
606 
607  Info<< "Initial nrmOrientation after comparison to NURBS u="
608  << u << ",v=" << v << " nrm: " << nrmOrientation_
609  << endl;
610 }
611 
612 
613 void NURBS3DSurface::flipNrmOrientation()
614 {
615  if (nrmOrientation_ == ALIGNED)
616  {
617  nrmOrientation_ = OPPOSED;
618  }
619  else
620  {
621  nrmOrientation_ = ALIGNED;
622  }
623 }
624 
625 
626 void NURBS3DSurface::setCPs(const List<vector>& CPs)
627 {
628  CPs_ = CPs;
629 }
630 
631 
632 void NURBS3DSurface::setWeights(const scalarList& weights)
633 {
634  weights_ = weights;
635 }
636 
637 
638 void NURBS3DSurface::setName(const word& name)
639 {
640  name_ = name;
641 }
642 
643 
644 void NURBS3DSurface::buildSurface()
645 {
646  const label uDegree(uBasis_.degree());
647  const label vDegree(vBasis_.degree());
648  const label uNCPs(uBasis_.nCPs());
649  const label vNCPs(vBasis_.nCPs());
650 /*
651  Info<< "\nuDegree: " << uDegree << "\nvDegree: " << vDegree
652  << "\nnUPts: " << nUPts_ << "\nnVPts: " << nVPts_
653  << "\nuNCPs: " << uNCPs << "\nvNCPs: " << vNCPs
654  << "\nNURBSSurface:\nCPs: " << CPs
655  << endl;
656 //*/
657  vectorField& field = *this;
658  field = vector::zero;
659 
660  for (label uI = 0; uI<nUPts_; uI++)
661  {
662  for (label vI = 0; vI<nVPts_; vI++)
663  {
664  const label ptI(uI*nVPts_ + vI);
665  const scalar& u(u_[ptI]);
666  const scalar& v(v_[ptI]);
667  scalar NMW(Zero);
668 
669  // Compute denominator.
670  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
671  {
672  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
673  {
674  const label CPI(vCPI*uNCPs + uCPI);
675 
676  NMW +=
677  uBasis_.basisValue(uCPI, uDegree, u)
678  * vBasis_.basisValue(vCPI, vDegree, v)
679  * weights_[CPI];
680  }
681  }
682 
683  // Compute the points.
684  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
685  {
686  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
687  {
688  const label CPI(vCPI*uNCPs + uCPI);
689 
690  this->operator[](ptI) +=
691  CPs_[CPI]
692  * uBasis_.basisValue(uCPI, uDegree, u)
693  * vBasis_.basisValue(vCPI, vDegree, v)
694  * weights_[CPI]/NMW;
695  }
696  }
697  }
698  }
699 }
700 
701 
702 
703 void NURBS3DSurface::invertU()
704 {
705  Info<< "Inverting NURBS surface " << name_ << " in u." << endl;
706 
707  const label uNCPs(uBasis_.nCPs());
708  const label vNCPs(vBasis_.nCPs());
709  List<vector> invertedCPs(CPs_.size(), Zero);
710  List<scalar> invertedWeights(CPs_.size(), Zero);
711 
712  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
713  {
714  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
715  {
716  const label CPI(vCPI*uNCPs + uCPI);
717  const label invUCPI(uNCPs-1-uCPI);
718  const label uInvCPI(vCPI*uNCPs + invUCPI);
719 
720  invertedCPs[CPI] = CPs_[uInvCPI];
721  invertedWeights[CPI] = weights_[uInvCPI];
722  }
723  }
724 
725  CPs_ = invertedCPs;
726  weights_ = invertedWeights;
727 
728  buildSurface();
729 }
730 
731 
732 void NURBS3DSurface::invertV()
733 {
734  Info<< "Inverting NURBS surface " << name_ << " in v." << endl;
735 
736  const label uNCPs(uBasis_.nCPs());
737  const label vNCPs(vBasis_.nCPs());
738  List<vector> invertedCPs(CPs_.size(), Zero);
739  List<scalar> invertedWeights(CPs_.size(), Zero);
740 
741  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
742  {
743  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
744  {
745  const label CPI(vCPI*uNCPs + uCPI);
746  const label invVCPI(vNCPs-1-vCPI);
747  const label vInvCPI(invVCPI*uNCPs + uCPI);
748 
749  invertedCPs[CPI] = CPs_[vInvCPI];
750  invertedWeights[CPI] = weights_[vInvCPI];
751  }
752  }
753 
754  CPs_ = invertedCPs;
755  weights_ = invertedWeights;
756 
757  buildSurface();
758 }
759 
760 
761 void NURBS3DSurface::invertUV()
762 {
763  Info<< "Inverting NURBS surface " << name_ << " in u and v." << endl;
764 
765  const label uNCPs(uBasis_.nCPs());
766  const label vNCPs(vBasis_.nCPs());
767  List<vector> invertedCPs(CPs_.size(), Zero);
768  List<scalar> invertedWeights(CPs_.size(), Zero);
769 
770  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
771  {
772  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
773  {
774  const label CPI(vCPI*uNCPs + uCPI);
775  const label invUCPI(uNCPs - 1 - uCPI);
776  const label invVCPI(vNCPs - 1 - vCPI);
777  const label uvInvCPI(invVCPI*uNCPs + invUCPI);
778 
779  invertedCPs[CPI] = CPs_[uvInvCPI];
780  invertedWeights[CPI] = weights_[uvInvCPI];
781  }
782  }
783 
784  CPs_ = invertedCPs;
785  weights_ = invertedWeights;
786 
787  buildSurface();
788 }
789 
790 
791 void NURBS3DSurface::makeEquidistant
792 (
793  const label lenAcc,
794  const label maxIter,
795  const label spacingCorrInterval,
796  const scalar tolerance
797 )
798 {
799 /*
800  Info<< "Making points equidistant is physical space on surface "
801  << name_
802  << endl;
803 //*/
804  // Equidistant spacing in u along v isoLines.
805  for (label vI = 0; vI<nVPts_; vI++)
806  {
807  scalarList UI(nUPts_, Zero);
808  const scalar VHeld(v_[vI]);
809  labelList uAddressing(nUPts_, -1);
810 
811  // Set the point uAddressing to re-assign correct u_ values later.
812  forAll(uAddressing, uI)
813  {
814  const label ptI(uI*nVPts_ + vI);
815  uAddressing[uI] = ptI;
816  }
817  // Set equidistant u values.
818  setEquidistantR
819  (
820  UI,
821  VHeld,
822  PARAMU,
823  lenAcc,
824  maxIter,
825  spacingCorrInterval,
826  tolerance
827  );
828 
829  // Re-assign new equidistant u values.
830  forAll(UI, uI)
831  {
832  const label& uAddress(uAddressing[uI]);
833  u_[uAddress] = UI[uI];
834  }
835  }
836 
837  // Equidistant spacing in v along u isoLines.
838  for (label uI = 0; uI<nUPts_; uI++)
839  {
840  scalarList VI(nVPts_, Zero);
841  const scalar UHeld(u_[uI*nVPts_]);
842  labelList vAddressing(nUPts_, -1);
843 
844  // Set the point vAddressing to re-assign correct u_ values later.
845  forAll(vAddressing, vI)
846  {
847  const label ptI(uI*nVPts_ + vI);
848  vAddressing[vI] = ptI;
849  }
850 
851  // Set equidistant u values.
852  setEquidistantR
853  (
854  VI,
855  UHeld,
856  PARAMV,
857  lenAcc,
858  maxIter,
859  spacingCorrInterval,
860  tolerance
861  );
862 
863  // Re-assign new equidistant u values.
864  forAll(VI, vI)
865  {
866  const label& vAddress(vAddressing[vI]);
867  v_[vAddress] = VI[vI];
868  }
869  }
870 
871  buildSurface();
872 }
873 
874 
875 // * * * * * * * * * * * * * Point Calc Functions * * * * * * * * * * * * * //
876 
877 vector NURBS3DSurface::surfacePoint
878 (
879  const scalar& u,
880  const scalar& v
881 )
882 {
883  const label uDegree(uBasis_.degree());
884  const label vDegree(vBasis_.degree());
885  const label uNCPs(uBasis_.nCPs());
886  const label vNCPs(vBasis_.nCPs());
887  scalar NMW(Zero);
888 
889  // Compute denominator.
890  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
891  {
892  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
893  {
894  const label CPI(vCPI*uNCPs + uCPI);
895 
896  NMW +=
897  uBasis_.basisValue(uCPI, uDegree, u)
898  * vBasis_.basisValue(vCPI, vDegree, v)
899  * weights_[CPI];
900  }
901  }
902 
903  // Compute the points.
904  vector point(Zero);
905 
906  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
907  {
908  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
909  {
910  const label CPI(vCPI*uNCPs + uCPI);
911 
912  point +=
913  CPs_[CPI]
914  * uBasis_.basisValue(uCPI, uDegree, u)
915  * vBasis_.basisValue(vCPI, vDegree, v)
916  * weights_[CPI]/NMW;
917  }
918  }
919 
920  return point;
921 }
922 
923 
924 scalarList NURBS3DSurface::findClosestSurfacePoint
925 (
926  const vector& targetPoint,
927  const label maxIter,
928  const scalar tolerance
929 )
930 {
931  // Loop through surface points to find the closest one for initialization.
932  const vectorField& surface(*this);
933  scalar dist(GREAT);
934  label closePtI(-1);
935 
936  forAll(surface, ptI)
937  {
938  const scalar distLoc(mag(targetPoint-surface[ptI]));
939 
940  if (distLoc < dist)
941  {
942  dist = distLoc;
943  closePtI = ptI;
944  }
945  }
946 
947  label iter(0);
948  scalar u(u_[closePtI]);
949  scalar v(v_[closePtI]);
950  vector xuv(surfacePoint(u, v));
951  scalar res(GREAT);
952  scalar resOld(GREAT);
953  scalar resDeriv(GREAT);
954  label nBoundsU(0);
955  label nBoundsV(0);
956 
957  do
958  {
959  /*
960  const vector dxdu(surfaceDerivativeU(u, v));
961  const vector dxdv(surfaceDerivativeV(u, v));
962  const vector d2xdu2(surfaceDerivativeUU(u, v));
963  const vector d2xdv2(surfaceDerivativeVV(u, v));
964  const scalar uLHS((dxdu&dxdu) + ((xuv-targetPoint) & d2xdu2));
965  const scalar uRHS(-((xuv-targetPoint) & dxdu));
966  const scalar vLHS((dxdv&dxdv) + ((xuv-targetPoint) & d2xdv2));
967  const scalar vRHS(-((xuv-targetPoint) & dxdv));
968 
969  // Update parametric coordinate and compute new point and
970  // bound param coordinates if needed.
971  // Compute residual.
972  u += uRHS/(uLHS+SMALL);
973  v += vRHS/(vLHS+SMALL);
974 
975  bound(u, v);
976 
977  xuv = surfacePoint(u, v);
978  res = mag((xuv-targetPoint) & surfaceDerivativeU(u, v))
979  + mag((xuv-targetPoint) & surfaceDerivativeV(u, v));
980  */
981 
982  const vector dxdu(surfaceDerivativeU(u, v));
983  const vector dxdv(surfaceDerivativeV(u, v));
984  const vector d2xdu2(surfaceDerivativeUU(u, v));
985  const vector d2xdv2(surfaceDerivativeVV(u, v));
986  const vector d2xduv(surfaceDerivativeUV(u, v));
987  const scalar a((dxdu&dxdu) + ((xuv-targetPoint) & d2xdu2));
988  const scalar b((dxdu&dxdv) + ((xuv-targetPoint) & d2xduv));
989  const scalar c=b;
990  const scalar d((dxdv&dxdv) + ((xuv-targetPoint) & d2xdv2));
991  const scalar invDenom = 1./(a*d-b*c);
992 
993  const scalar uRHS(-((xuv-targetPoint) & dxdu));
994  const scalar vRHS(-((xuv-targetPoint) & dxdv));
995 
996  // Update parametric coordinate and compute new point and
997  // bound param coordinates if needed.
998  // Compute residual.
999  u += ( d*uRHS-b*vRHS)*invDenom;
1000  v += (-c*uRHS+a*vRHS)*invDenom;
1001 
1002  if (boundDirection(u))
1003  {
1004  nBoundsU++;
1005  }
1006  if (boundDirection(v))
1007  {
1008  nBoundsV++;
1009  }
1010 
1011  xuv = surfacePoint(u, v);
1012  // If manual assignment in u is required, deal only with the v eqn
1013  if (nBoundsU >= 5)
1014  {
1015  res = mag((xuv-targetPoint) & surfaceDerivativeV(u, v));
1016  resDeriv = mag(res-resOld)/resOld;
1017  resOld = res;
1018 // Info<< targetPoint << " " << res << endl;
1019  }
1020  // If manual assignment in v is required, deal only with the u eqn
1021  else if (nBoundsV >= 5)
1022  {
1023  res = mag((xuv-targetPoint) & surfaceDerivativeU(u, v));
1024  resDeriv = mag(res-resOld)/resOld;
1025  resOld = res;
1026 // Info<< targetPoint << " " << res << endl;
1027  }
1028  else if (nBoundsU <= 5 && nBoundsV <= 5)
1029  {
1030  res = mag((xuv-targetPoint) & surfaceDerivativeU(u, v))
1031  + mag((xuv-targetPoint) & surfaceDerivativeV(u, v));
1032  resDeriv = mag(res-resOld)/resOld;
1033  resOld = res;
1034  }
1035  else
1036  {
1037  WarningInFunction
1038  << "More than 5 bounds in both the u and v directions!"
1039  << "Something seems weird" << nBoundsU << " " << nBoundsV
1040  << endl;
1041  res = mag((xuv-targetPoint) & surfaceDerivativeU(u, v))
1042  + mag((xuv-targetPoint) & surfaceDerivativeV(u, v));
1043  resDeriv = mag(res-resOld)/resOld;
1044  resOld = res;
1045  }
1046  }
1047  while ((iter++ < maxIter) && (res > tolerance));
1048 
1049  scalarList closestParameters(2, u);
1050  closestParameters[1] = v;
1051  // warning if method did not reach threshold
1052  if (iter > maxIter)
1053  {
1054  WarningInFunction
1055  << "Finding surface point closest to " << targetPoint
1056  << " for surface " << name_ << " failed \n"
1057  << " Number of bounding operations in u,v "
1058  << nBoundsU << " " << nBoundsV << endl
1059  << " Residual value and derivative " << res << " " << resDeriv
1060  << endl << endl;
1061  closestParameters = -1;
1062  }
1063 
1064  return closestParameters;
1065 }
1066 
1067 
1068 tmp<vector2DField> NURBS3DSurface::findClosestSurfacePoint
1069 (
1070  const vectorField& targetPoints,
1071  const label maxIter,
1072  const scalar tolerance
1073 )
1074 {
1075  auto tparamCoors = tmp<vector2DField>::New(targetPoints.size(), Zero);
1076  auto& paramCoors = tparamCoors.ref();
1077 
1078  const vectorField& surface(*this);
1079  label nBoundedPoints(0);
1080  scalar maxResidual(0);
1081  scalar maxResidualDeriv(0);
1082  forAll(paramCoors, pI)
1083  {
1084  const vector& targetPoint(targetPoints[pI]);
1085 
1086  // Loop through surface points to find the closest one for
1087  // initialization. The initialization could possibly be done with a
1088  // geodesical Laplace. Potentially faster?
1089  scalar dist(GREAT);
1090  label closePtI(-1);
1091 
1092  forAll(surface, ptI)
1093  {
1094  const scalar distLoc(mag(targetPoint - surface[ptI]));
1095 
1096  if (distLoc < dist)
1097  {
1098  dist = distLoc;
1099  closePtI = ptI;
1100  }
1101  }
1102 
1103  label iter(0);
1104  scalar u(u_[closePtI]);
1105  scalar v(v_[closePtI]);
1106  vector xuv(surfacePoint(u, v));
1107  scalar res(GREAT);
1108  scalar resOld(GREAT);
1109  scalar resDeriv(GREAT);
1110  label nBoundsU(0);
1111  label nBoundsV(0);
1112 
1113  do
1114  {
1115  const vector dxdu(surfaceDerivativeU(u, v));
1116  const vector dxdv(surfaceDerivativeV(u, v));
1117  const vector d2xdu2(surfaceDerivativeUU(u, v));
1118  const vector d2xdv2(surfaceDerivativeVV(u, v));
1119  const vector d2xduv(surfaceDerivativeUV(u, v));
1120  const scalar a((dxdu&dxdu) + ((xuv-targetPoint) & d2xdu2));
1121  const scalar b((dxdu&dxdv) + ((xuv-targetPoint) & d2xduv));
1122  const scalar c=b;
1123  const scalar d((dxdv&dxdv) + ((xuv-targetPoint) & d2xdv2));
1124  const scalar invDenom = 1./(a*d-b*c);
1125 
1126  const scalar uRHS(-((xuv-targetPoint) & dxdu));
1127  const scalar vRHS(-((xuv-targetPoint) & dxdv));
1128 
1129  // Update parametric coordinate and compute new point and
1130  // bound param coordinates if needed.
1131  // Compute residual.
1132  u += ( d*uRHS-b*vRHS)*invDenom;
1133  v += (-c*uRHS+a*vRHS)*invDenom;
1134 
1135  if (boundDirection(u))
1136  {
1137  nBoundsU++;
1138  }
1139  if (boundDirection(v))
1140  {
1141  nBoundsV++;
1142  }
1143 
1144  xuv = surfacePoint(u, v);
1145  // If manual assignment in u is required, deal only with the v eqn
1146  if (nBoundsU >= 5)
1147  {
1148  res = mag((xuv-targetPoint) & surfaceDerivativeV(u, v));
1149  resDeriv = mag(res-resOld)/resOld;
1150  resOld = res;
1151  }
1152  // If manual assignment in b is required, deal only with the u eqn
1153  else if (nBoundsV >= 5)
1154  {
1155  res = mag((xuv-targetPoint) & surfaceDerivativeU(u, v));
1156  resDeriv = mag(res-resOld)/resOld;
1157  resOld = res;
1158  }
1159  else if (nBoundsU<=5 && nBoundsV<=5)
1160  {
1161  res = mag((xuv-targetPoint) & surfaceDerivativeU(u, v))
1162  + mag((xuv-targetPoint) & surfaceDerivativeV(u, v));
1163  resDeriv = mag(res-resOld)/resOld;
1164  resOld = res;
1165  }
1166  else
1167  {
1168  WarningInFunction
1169  << "More than 5 bounding operations in both the u and v directions!"
1170  << "u direction " << nBoundsU << endl
1171  << "v direction " << nBoundsV << endl
1172  << endl;
1173  res =
1174  mag((xuv-targetPoint) & surfaceDerivativeU(u, v))
1175  + mag((xuv-targetPoint) & surfaceDerivativeV(u, v));
1176  resDeriv = mag(res-resOld)/resOld;
1177  resOld = res;
1178  }
1179  }
1180  while ((iter++ < maxIter) && (res > tolerance));
1181 
1182  // warning if method did not reach threshold
1183  if (iter > maxIter)
1184  {
1185  nBoundedPoints++;
1186  maxResidual = max(maxResidual, res);
1187  maxResidualDeriv = max(maxResidualDeriv, resDeriv);
1188  }
1189 
1190  paramCoors[pI].x() = u;
1191  paramCoors[pI].y() = v;
1192  }
1193  reduce(nBoundedPoints, sumOp<label>());
1194  reduce(maxResidual, maxOp<scalar>());
1195  reduce(maxResidualDeriv, maxOp<scalar>());
1196  Info<< "Points that couldn't reach the residual limit:: "
1197  << nBoundedPoints << endl
1198  << "Max residual after reaching iterations limit:: "
1199  << maxResidual << endl
1200  << "Max residual derivative after reaching iterations limit:: "
1201  << maxResidualDeriv << endl
1202  << endl;
1203 
1204  return tparamCoors;
1205 }
1206 
1207 
1208 scalarList NURBS3DSurface::findClosestSurfacePoint
1209 (
1210  const vector& targetPoint,
1211  const scalar& uInitGuess,
1212  const scalar& vInitGuess,
1213  const label maxIter,
1214  const scalar tolerance
1215 )
1216 {
1217  // Loop through surface points to find the closest one for initialization.
1218  label iter(0);
1219  scalar u(uInitGuess);
1220  scalar v(vInitGuess);
1221  vector xuv(surfacePoint(u, v));
1222  scalar res(GREAT);
1223 
1224  do
1225  {
1226  const vector dxdu(surfaceDerivativeU(u, v));
1227  const vector dxdv(surfaceDerivativeV(u, v));
1228  const vector d2xdu2(surfaceDerivativeUU(u, v));
1229  const vector d2xdv2(surfaceDerivativeVV(u, v));
1230  const scalar uLHS((dxdu&dxdu) + ((xuv-targetPoint) & d2xdu2));
1231  const scalar uRHS(-((xuv-targetPoint) & dxdu));
1232  const scalar vLHS((dxdv&dxdv) + ((xuv-targetPoint) & d2xdv2));
1233  const scalar vRHS(-((xuv-targetPoint) & dxdv));
1234 
1235  // Update parametric coordinate and compute new point and
1236  // bound param coordinates if needed.
1237  // Compute residual.
1238  u += uRHS/(uLHS+SMALL);
1239  v += vRHS/(vLHS+SMALL);
1240 
1241  bound(u, v);
1242 
1243  xuv = surfacePoint(u, v);
1244  res =
1245  mag((xuv-targetPoint) & surfaceDerivativeU(u, v))
1246  + mag((xuv-targetPoint) & surfaceDerivativeV(u, v));
1247  }
1248  while ((iter++ < maxIter) && (res > tolerance));
1249 
1250  // warning if method did not reach threshold
1251  if (iter > maxIter)
1252  {
1253  WarningInFunction
1254  << "Finding surface point closest to " << targetPoint << " failed."
1255  << endl;
1256  }
1257 
1258  scalarList closestParameters(2, u);
1259  closestParameters[1] = v;
1260 
1261  return closestParameters;
1262 }
1263 
1264 
1265 const vector NURBS3DSurface::nrm(scalar u, scalar v)
1266 {
1267  vector surfaceNrm(Zero);
1268 
1269  if (nrmOrientation_ == ALIGNED)
1270  {
1271  surfaceNrm = surfaceDerivativeU(u, v) ^ surfaceDerivativeV(u, v);
1272  }
1273  else
1274  {
1275  surfaceNrm = surfaceDerivativeV(u, v) ^ surfaceDerivativeU(u, v);
1276  }
1277 
1278  surfaceNrm /= mag(surfaceNrm);
1279 
1280  return surfaceNrm;
1281 }
1282 
1283 
1284 List<scalarList> NURBS3DSurface::genEquidistant
1285 (
1286  const label nUPts,
1287  const label nVPts,
1288  const label lenAcc,
1289  const label maxIter,
1290  const label spacingCorrInterval,
1291  const scalar tolerance
1292 )
1293 {
1294 /*
1295  Info<< "Generating points equidistant in physical space on surface "
1296  << name_
1297  << endl;
1298 //*/
1299  // Preset U and V with uniform values.
1300  List<scalarList> UV(NPARAMS, scalarList(0));
1301 
1302  scalarList& U(UV[PARAMU]);
1303  scalarList& V(UV[PARAMV]);
1304 
1305  setUniformUV(U, V, nUPts, nVPts);
1306 
1307  // Equidistant spacing in u along v isoLines.
1308  for (label vI = 0; vI<nVPts; vI++)
1309  {
1310  scalarList UI(nUPts, Zero);
1311  const scalar VHeld(V[vI]);
1312  labelList uAddressing(nUPts, -1);
1313 
1314  // Set the point uAddressing to re-assign correct U values later.
1315  forAll(uAddressing, uI)
1316  {
1317  const label ptI(uI*nVPts + vI);
1318  uAddressing[uI] = ptI;
1319  }
1320  // Set equidistant u values.
1321  setEquidistantR
1322  (
1323  UI,
1324  VHeld,
1325  PARAMU,
1326  lenAcc,
1327  maxIter,
1328  spacingCorrInterval,
1329  tolerance
1330  );
1331 
1332  // Re-assign new equidistant u values.
1333  forAll(UI, uI)
1334  {
1335  const label& uAddress(uAddressing[uI]);
1336  U[uAddress] = UI[uI];
1337  }
1338  }
1339 
1340  // Equidistant spacing in v along u isoLines.
1341  for (label uI = 0; uI<nUPts; uI++)
1342  {
1343  scalarList VI(nVPts, Zero);
1344  const scalar UHeld(U[uI*nVPts]);
1345  labelList vAddressing(nUPts, -1);
1346 
1347  // Set the point vAddressing to re-assign correct V values later.
1348  forAll(vAddressing, vI)
1349  {
1350  const label ptI(uI*nVPts + vI);
1351  vAddressing[vI] = ptI;
1352  }
1353 
1354  // Set equidistant u values.
1355  setEquidistantR
1356  (
1357  VI,
1358  UHeld,
1359  PARAMV,
1360  lenAcc,
1361  maxIter,
1362  spacingCorrInterval,
1363  tolerance
1364  );
1365 
1366  // Re-assign new equidistant u values.
1367  forAll(VI, vI)
1368  {
1369  const label& vAddress(vAddressing[vI]);
1370  V[vAddress] = VI[vI];
1371  }
1372  }
1373 
1374  return UV;
1375 }
1376 
1377 
1378 // * * * * * * * * * * * * * * Location Functions * * * * * * * * * * * * * //
1379 
1380 bool NURBS3DSurface::checkRangeU
1381 (
1382  const scalar u,
1383  const label CPI,
1384  const label uDegree
1385 ) const
1386 {
1387  const label uCPI(CPsUCPIs_[CPI]);
1388 
1389  return uBasis_.checkRange(u, uCPI, uDegree);
1390 }
1391 
1392 
1393 bool NURBS3DSurface::checkRangeU
1394 (
1395  const scalar u,
1396  const label CPI
1397 ) const
1398 {
1399  const label uDegree(uBasis_.degree());
1400 
1401  return checkRangeU(u, CPI, uDegree);
1402 }
1403 
1404 
1405 bool NURBS3DSurface::checkRangeV
1406 (
1407  const scalar v,
1408  const label CPI,
1409  const label vDegree
1410 ) const
1411 {
1412  const label vCPI(CPsVCPIs_[CPI]);
1413 
1414  return vBasis_.checkRange(v, vCPI, vDegree);
1415 }
1416 
1417 
1418 bool NURBS3DSurface::checkRangeV
1419 (
1420  const scalar v,
1421  const label CPI
1422 ) const
1423 {
1424  const label vDegree(vBasis_.degree());
1425 
1426  return checkRangeV(v, CPI, vDegree);
1427 }
1428 
1429 
1430 bool NURBS3DSurface::checkRangeUV
1431 (
1432  const scalar v,
1433  const scalar u,
1434  const label CPI,
1435  const label uDegree,
1436  const label vDegree
1437 ) const
1438 {
1439  if (checkRangeU(u, CPI, uDegree) && checkRangeV(v, CPI, vDegree))
1440  {
1441  return true;
1442  }
1443 
1444  return false;
1445 }
1446 
1447 
1448 bool NURBS3DSurface::checkRangeUV
1449 (
1450  const scalar v,
1451  const scalar u,
1452  const label CPI
1453 ) const
1454 {
1455  const label uDegree(uBasis_.degree());
1456  const label vDegree(vBasis_.degree());
1457 
1458  return checkRangeUV(u, v, CPI, uDegree, vDegree);
1459 }
1460 
1461 
1462 scalar NURBS3DSurface::lengthU
1463 (
1464  const label vIConst,
1465  const label uIStart,
1466  const label uIEnd
1467 ) const
1468 {
1469  // Compute derivatives wrt u for the given u interval.
1470  const label uLenSize(uIEnd - uIStart + 1);
1471  vectorField dxdu(uLenSize, Zero);
1472  scalar uLength(Zero);
1473 
1474  forAll(dxdu, uI)
1475  {
1476  const label ptI((uIStart+uI)*nVPts_ + vIConst);
1477  const label& u(u_[ptI]);
1478  const label& v(v_[ptI]);
1479 
1480  dxdu[uI] = surfaceDerivativeU(u, v);
1481  }
1482 
1483  // Integrate.
1484  for(label uI = 0; uI<(uLenSize - 1); uI++)
1485  {
1486  const label ptI((uIStart+uI)*nVPts_ + vIConst);
1487 
1488  uLength +=
1489  0.5*(mag(dxdu[uI + 1]) + mag(dxdu[uI]))*(u_[ptI + 1]-u_[ptI]);
1490  }
1491 
1492  return uLength;
1493 }
1494 
1495 
1496 scalar NURBS3DSurface::lengthU
1497 (
1498  const scalar vConst,
1499  const scalar uStart,
1500  const scalar uEnd,
1501  const label nPts
1502 ) const
1503 {
1504  // Compute derivatives wrt u for the given u interval.
1505  vectorField dxdu(nPts, Zero);
1506  scalarField localU(nPts, Zero);
1507  scalar uLength(Zero);
1508 
1509  forAll(localU, uI)
1510  {
1511  scalar& uLocal(localU[uI]);
1512  uLocal = uStart + scalar(uI)/scalar(nPts-1)*(uEnd-uStart);
1513  dxdu[uI] = surfaceDerivativeU(uLocal, vConst);
1514  }
1515 
1516  // Integrate.
1517  for(label uI = 0; uI<(nPts - 1); uI++)
1518  {
1519  uLength +=
1520  0.5*(mag(dxdu[uI + 1]) + mag(dxdu[uI]))*(localU[uI + 1]-localU[uI]);
1521  }
1522 
1523  return uLength;
1524 }
1525 
1526 
1527 scalar NURBS3DSurface::lengthU(const label vIConst) const
1528 {
1529  return lengthU(vIConst, 0, (nUPts_ - 1));
1530 }
1531 
1532 
1533 scalar NURBS3DSurface::lengthU(const scalar vConst) const
1534 {
1535  return lengthU(vConst, scalar(0), scalar(1), 100);
1536 }
1537 
1538 
1539 scalar NURBS3DSurface::lengthV
1540 (
1541  const label uIConst,
1542  const label vIStart,
1543  const label vIEnd
1544 ) const
1545 {
1546  // Compute derivatives wrt v for the given v interval.
1547  const label vLenSize(vIEnd - vIStart + 1);
1548  vectorField dxdv(vLenSize, Zero);
1549  scalar vLength(Zero);
1550 
1551  forAll(dxdv, vI)
1552  {
1553  const label ptI((uIConst)*nVPts_ + (vIStart+vI));
1554  const label& u(u_[ptI]);
1555  const label& v(v_[ptI]);
1556 
1557  dxdv[vI] = surfaceDerivativeV(u, v);
1558  }
1559 
1560  // Integrate.
1561  for(label vI = 0; vI<(vLenSize - 1); vI++)
1562  {
1563  const label ptI((uIConst)*nVPts_ + (vIStart + vI));
1564 
1565  vLength +=
1566  0.5*(mag(dxdv[vI + 1]) + mag(dxdv[vI]))*(v_[ptI + 1] - v_[ptI]);
1567  }
1568 
1569  return vLength;
1570 }
1571 
1572 
1573 scalar NURBS3DSurface::lengthV
1574 (
1575  const scalar uConst,
1576  const scalar vStart,
1577  const scalar vEnd,
1578  const label nPts
1579 ) const
1580 {
1581  // Compute derivatives wrt v for the given v interval.
1582  vectorField dxdv(nPts, Zero);
1583  scalarField localV(nPts, Zero);
1584  scalar vLength(Zero);
1585 
1586  forAll(localV, vI)
1587  {
1588  scalar& vLocal(localV[vI]);
1589  vLocal = vStart + scalar(vI)/scalar(nPts - 1)*(vEnd - vStart);
1590  dxdv[vI] = surfaceDerivativeV(uConst, vLocal);
1591  }
1592 
1593  // Integrate.
1594  for(label vI = 0; vI<(nPts - 1); vI++)
1595  {
1596  vLength +=
1597  0.5*(mag(dxdv[vI + 1]) + mag(dxdv[vI]))*(localV[vI + 1]-localV[vI]);
1598  }
1599 
1600  return vLength;
1601 }
1602 
1603 
1604 scalar NURBS3DSurface::lengthV(const label uIConst) const
1605 {
1606  return lengthV(uIConst, 0, (nVPts_ - 1));
1607 }
1608 
1609 
1610 scalar NURBS3DSurface::lengthV(const scalar uConst) const
1611 {
1612  return lengthV(uConst, scalar(0), scalar(1), 100);
1613 }
1614 
1615 
1616 // * * * * * * * * * * * * * Derivative Functions * * * * * * * * * * * * * //
1617 
1618 vector NURBS3DSurface::surfaceDerivativeU
1619 (
1620  const scalar uIn,
1621  const scalar vIn
1622 ) const
1623 {
1624  const label uDegree(uBasis_.degree());
1625  const label vDegree(vBasis_.degree());
1626  const label uNCPs(uBasis_.nCPs());
1627  const label vNCPs(vBasis_.nCPs());
1628  vector NMWP(Zero);
1629  vector dNduMWP(Zero);
1630  scalar NMW(Zero);
1631  scalar dNduMW(Zero);
1632 
1633  scalar u = uIn;
1634  scalar v = vIn;
1635  bound(u, v);
1636 
1637  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
1638  {
1639  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
1640  {
1641  const label CPI(vCPI*uNCPs + uCPI);
1642  const scalar uBasisValue(uBasis_.basisValue(uCPI, uDegree, u));
1643  const scalar vBasisValue(vBasis_.basisValue(vCPI, vDegree, v));
1644  const scalar uBasisDeriv
1645  (uBasis_.basisDerivativeU(uCPI, uDegree, u));
1646 
1647  NMWP += uBasisValue * vBasisValue * weights_[CPI] * CPs_[CPI];
1648  dNduMWP += uBasisDeriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1649  NMW += uBasisValue * vBasisValue * weights_[CPI];
1650  dNduMW += uBasisDeriv * vBasisValue * weights_[CPI];
1651  }
1652  }
1653 
1654  const vector uDerivative((dNduMWP - dNduMW*NMWP/(NMW+SMALL))/(NMW+SMALL));
1655 
1656  return uDerivative;
1657 }
1658 
1659 
1660 vector NURBS3DSurface::surfaceDerivativeV
1661 (
1662  const scalar uIn,
1663  const scalar vIn
1664 ) const
1665 {
1666  const label uDegree(uBasis_.degree());
1667  const label vDegree(vBasis_.degree());
1668  const label uNCPs(uBasis_.nCPs());
1669  const label vNCPs(vBasis_.nCPs());
1670  vector NMWP(Zero);
1671  vector dMdvNWP(Zero);
1672  scalar NMW(Zero);
1673  scalar dMdvNW(Zero);
1674 
1675  scalar u = uIn;
1676  scalar v = vIn;
1677  bound(u, v);
1678 
1679  for (label vCPI = 0; vCPI<vNCPs; vCPI++)
1680  {
1681  for (label uCPI = 0; uCPI<uNCPs; uCPI++)
1682  {
1683  const label CPI(vCPI*uNCPs + uCPI);
1684  const scalar uBasisValue(uBasis_.basisValue(uCPI, uDegree, u));
1685  const scalar vBasisValue(vBasis_.basisValue(vCPI, vDegree, v));
1686  const scalar vBasisDeriv
1687  (vBasis_.basisDerivativeU(vCPI, vDegree, v));
1688 
1689  NMWP += uBasisValue * vBasisValue * weights_[CPI] * CPs_[CPI];
1690  dMdvNWP += vBasisDeriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1691  NMW += uBasisValue * vBasisValue * weights_[CPI];
1692  dMdvNW += vBasisDeriv * uBasisValue * weights_[CPI];
1693  }
1694  }
1695 
1696  const vector vDerivative((dMdvNWP - dMdvNW*NMWP/(NMW+SMALL))/(NMW+SMALL));
1697 
1698  return vDerivative;
1699 }
1700 
1701 
1702 vector NURBS3DSurface::surfaceDerivativeUV
1703 (
1704  const scalar uIn,
1705  const scalar vIn
1706 ) const
1707 {
1708  const label uDegree(uBasis_.degree());
1709  const label vDegree(vBasis_.degree());
1710  const label uNCPs(uBasis_.nCPs());
1711  const label vNCPs(vBasis_.nCPs());
1712  vector NMWP(Zero);
1713  vector dNduMWP(Zero);
1714  vector dMdvNWP(Zero);
1715  vector dNMduvWP(Zero);
1716  scalar NMW(Zero);
1717  scalar dNduMW(Zero);
1718  scalar dMdvNW(Zero);
1719  scalar dNMduvW(Zero);
1720 
1721  scalar u = uIn;
1722  scalar v = vIn;
1723  bound(u, v);
1724 
1725  for (label vCPI=0; vCPI<vNCPs; vCPI++)
1726  {
1727  for (label uCPI=0; uCPI<uNCPs; uCPI++)
1728  {
1729  const label CPI(vCPI*uNCPs + uCPI);
1730  const scalar uBasisValue(uBasis_.basisValue(uCPI, uDegree, u));
1731  const scalar vBasisValue(vBasis_.basisValue(vCPI, vDegree, v));
1732  const scalar uBasisDeriv
1733  (uBasis_.basisDerivativeU(uCPI, uDegree, u));
1734  const scalar vBasisDeriv
1735  (vBasis_.basisDerivativeU(vCPI, vDegree, v));
1736  //Info<< "vCPI=" << vCPI << ",uCPI=" << uCPI
1737  // << " N=" << uBasisValue << ",N'=" << uBasisDeriv
1738  // << " M=" << vBasisValue << ",M'=" << vBasisDeriv
1739  // << endl;
1740  NMWP += vBasisValue * uBasisValue * weights_[CPI] * CPs_[CPI];
1741  dNduMWP += uBasisDeriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1742  dMdvNWP += vBasisDeriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1743  dNMduvWP += uBasisDeriv * vBasisDeriv * weights_[CPI] * CPs_[CPI];
1744  NMW += vBasisValue * uBasisValue * weights_[CPI];
1745  dNduMW += uBasisDeriv * vBasisValue * weights_[CPI];
1746  dMdvNW += vBasisDeriv * uBasisValue * weights_[CPI];
1747  dNMduvW += uBasisDeriv * vBasisDeriv * weights_[CPI];
1748  }
1749  }
1750 
1751  const vector uvDerivative
1752  (
1753  (
1754  dNMduvWP
1755  - (dNMduvW*NMWP + dMdvNW*dNduMWP + dNduMW*dMdvNWP)/(NMW+SMALL)
1756  + scalar(2)*dNduMW*dMdvNW*NMWP/(NMW+SMALL)/(NMW+SMALL)
1757  ) / (NMW+SMALL)
1758  );
1759 
1760  return uvDerivative;
1761 }
1762 
1763 
1764 vector NURBS3DSurface::surfaceDerivativeUU
1765 (
1766  const scalar uIn,
1767  const scalar vIn
1768 ) const
1769 {
1770  const label uDegree(uBasis_.degree());
1771  const label vDegree(vBasis_.degree());
1772  const label uNCPs(uBasis_.nCPs());
1773  const label vNCPs(vBasis_.nCPs());
1774  vector NMWP(Zero);
1775  vector dNduMWP(Zero);
1776  vector d2Ndu2MWP(Zero);
1777  scalar NMW(Zero);
1778  scalar dNduMW(Zero);
1779  scalar d2Ndu2MW(Zero);
1780 
1781  scalar u = uIn;
1782  scalar v = vIn;
1783  bound(u, v);
1784 
1785  for (label vCPI=0; vCPI<vNCPs; vCPI++)
1786  {
1787  for (label uCPI=0; uCPI<uNCPs; uCPI++)
1788  {
1789  const label CPI(vCPI*uNCPs + uCPI);
1790  const scalar uBasisValue(uBasis_.basisValue(uCPI, uDegree, u));
1791  const scalar vBasisValue(vBasis_.basisValue(vCPI, vDegree, v));
1792  const scalar uBasisDeriv
1793  (uBasis_.basisDerivativeU(uCPI, uDegree, u));
1794  const scalar uBasis2Deriv
1795  (uBasis_.basisDerivativeUU(uCPI, uDegree, u));
1796 
1797  NMWP += uBasisValue * vBasisValue * weights_[CPI] * CPs_[CPI];
1798  dNduMWP += uBasisDeriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1799  d2Ndu2MWP += uBasis2Deriv * vBasisValue * weights_[CPI] * CPs_[CPI];
1800  NMW += uBasisValue * vBasisValue * weights_[CPI];
1801  dNduMW += uBasisDeriv * vBasisValue * weights_[CPI];
1802  d2Ndu2MW += uBasis2Deriv * vBasisValue * weights_[CPI];
1803  }
1804  }
1805 
1806  const vector uuDerivative
1807  (
1808  (
1809  d2Ndu2MWP
1810  - scalar(2)*dNduMW*dNduMWP/(NMW+SMALL)
1811  - d2Ndu2MW*NMWP/(NMW+SMALL)
1812  + scalar(2)*dNduMW*dNduMW*NMWP/(NMW+SMALL)/(NMW+SMALL)
1813  ) / (NMW+SMALL)
1814  );
1815 
1816  return uuDerivative;
1817 }
1818 
1819 
1820 vector NURBS3DSurface::surfaceDerivativeVV
1821 (
1822  const scalar uIn,
1823  const scalar vIn
1824 ) const
1825 {
1826  const label uDegree(uBasis_.degree());
1827  const label vDegree(vBasis_.degree());
1828  const label uNCPs(uBasis_.nCPs());
1829  const label vNCPs(vBasis_.nCPs());
1830  vector NMWP(Zero);
1831  vector dMdvNWP(Zero);
1832  vector d2Mdv2NWP(Zero);
1833  scalar NMW(Zero);
1834  scalar dMdvNW(Zero);
1835  scalar d2Mdv2NW(Zero);
1836 
1837  scalar u = uIn;
1838  scalar v = vIn;
1839  bound(u, v);
1840 
1841  for (label vCPI=0; vCPI<vNCPs; vCPI++)
1842  {
1843  for (label uCPI=0; uCPI<uNCPs; uCPI++)
1844  {
1845  const label CPI(vCPI*uNCPs + uCPI);
1846  const scalar uBasisValue(uBasis_.basisValue(uCPI, uDegree, u));
1847  const scalar vBasisValue(vBasis_.basisValue(vCPI, vDegree, v));
1848  const scalar vBasisDeriv
1849  (vBasis_.basisDerivativeU(vCPI, vDegree, v));
1850  const scalar vBasis2Deriv
1851  (vBasis_.basisDerivativeUU(vCPI, vDegree, v));
1852 
1853  NMWP += vBasisValue * uBasisValue * weights_[CPI] * CPs_[CPI];
1854  dMdvNWP += vBasisDeriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1855  d2Mdv2NWP += vBasis2Deriv * uBasisValue * weights_[CPI] * CPs_[CPI];
1856  NMW += vBasisValue * uBasisValue * weights_[CPI];
1857  dMdvNW += vBasisDeriv * uBasisValue * weights_[CPI];
1858  d2Mdv2NW += vBasis2Deriv * uBasisValue * weights_[CPI];
1859  }
1860  }
1861 
1862  const vector vvDerivative
1863  (
1864  (
1865  d2Mdv2NWP
1866  - scalar(2)*dMdvNW*dMdvNWP/(NMW+SMALL)
1867  - d2Mdv2NW*NMWP/(NMW+SMALL)
1868  + scalar(2)*dMdvNW*dMdvNW*NMWP/(NMW+SMALL)/(NMW+SMALL)
1869  ) / (NMW+SMALL)
1870  );
1871 
1872  return vvDerivative;
1873 }
1874 
1875 
1876 scalar NURBS3DSurface::surfaceDerivativeCP
1877 (
1878  const scalar u,
1879  const scalar v,
1880  const label CPI
1881 ) const
1882 {
1883  //Info<< "u,v,cpI " << u << " " << v << " " << CPI << endl;
1884  // compute denominator.
1885  const label uDegree(uBasis_.degree());
1886  const label vDegree(vBasis_.degree());
1887  const label uNCPs(uBasis_.nCPs());
1888  const label vNCPs(vBasis_.nCPs());
1889  const label uCPI(CPsUCPIs_[CPI]);
1890  const label vCPI(CPsVCPIs_[CPI]);
1891  scalar NMW(Zero);
1892 
1893  for (label vCPJ=0; vCPJ<vNCPs; vCPJ++)
1894  {
1895  for (label uCPJ=0; uCPJ<uNCPs; uCPJ++)
1896  {
1897  const label CPJ(vCPJ*uNCPs + uCPJ);
1898  const scalar uBasisValue(uBasis_.basisValue(uCPJ, uDegree, u));
1899  const scalar vBasisValue(vBasis_.basisValue(vCPJ, vDegree, v));
1900 
1901  NMW += vBasisValue * uBasisValue * weights_[CPJ];
1902  }
1903  }
1904  //Info<< "denom " << NMW << endl;
1905 
1906  const scalar CPDerivative
1907  (
1908  uBasis_.basisValue(uCPI, uDegree, u)
1909  * vBasis_.basisValue(vCPI, vDegree, v)
1910  * weights_[CPI]
1911  / (NMW+SMALL)
1912  );
1913 
1914  return CPDerivative;
1915 }
1916 
1917 
1918 vector NURBS3DSurface::surfaceDerivativeW
1919 (
1920  const scalar u,
1921  const scalar v,
1922  const label CPI
1923 ) const
1924 {
1925  // Compute nominator, denominator.
1926  const label uDegree(uBasis_.degree());
1927  const label vDegree(vBasis_.degree());
1928  const label uNCPs(uBasis_.nCPs());
1929  const label vNCPs(vBasis_.nCPs());
1930  const label uCPI(CPsUCPIs_[CPI]);
1931  const label vCPI(CPsVCPIs_[CPI]);
1932  vector NMWP(Zero);
1933  scalar NMW(Zero);
1934 
1935  for (label vCPJ=0; vCPJ<vNCPs; vCPJ++)
1936  {
1937  for (label uCPJ=0; uCPJ<uNCPs; uCPJ++)
1938  {
1939  const label CPJ(vCPJ*uNCPs + uCPJ);
1940  const scalar uBasisValue(uBasis_.basisValue(uCPJ, uDegree, u));
1941  const scalar vBasisValue(vBasis_.basisValue(vCPJ, vDegree, v));
1942 
1943  NMWP += uBasisValue * vBasisValue * weights_[CPJ] * CPs_[CPJ];
1944  NMW += uBasisValue * vBasisValue * weights_[CPJ];
1945  }
1946  }
1947 
1948  // Compute derivative.
1949  const vector WDerivative
1950  (
1951  uBasis_.basisValue(uCPI, uDegree, u)
1952  * vBasis_.basisValue(vCPI, vDegree, v)
1953  * (CPs_[CPI] - (NMWP / (NMW+SMALL)))
1954  / (NMW+SMALL)
1955  );
1956 
1957  return WDerivative;
1958 }
1959 
1960 
1961 scalar NURBS3DSurface::lengthDerivativeU
1962 (
1963  const scalar vConst,
1964  const scalar uStart,
1965  const scalar uEnd,
1966  const label nPts
1967 ) const
1968 {
1969  // compute derivatives wrt u for the given u interval
1970  vectorField dxdu(nPts, Zero);
1971  vectorField d2xdu2(nPts, Zero);
1972  scalarField localU(nPts, Zero);
1973  scalar ulDerivative(Zero);
1974 
1975  forAll(localU, uI)
1976  {
1977  scalar& uLocal(localU[uI]);
1978  uLocal = uStart + scalar(uI)/scalar(nPts-1)*(uEnd-uStart);
1979  dxdu[uI] = surfaceDerivativeU(uLocal, vConst);
1980  d2xdu2[uI] = surfaceDerivativeUU(uLocal, vConst);
1981  }
1982 
1983  // Integrate.
1984  for(label uI=0; uI<(nPts-1); uI++)
1985  {
1986  ulDerivative +=
1987  0.5
1988  * (
1989  (dxdu[uI+1]&d2xdu2[uI+1])/(mag(dxdu[uI+1])+SMALL)
1990  + (dxdu[uI ]&d2xdu2[uI ])/(mag(dxdu[uI ])+SMALL)
1991  )
1992  * (localU[uI+1]-localU[uI]);
1993  }
1994 
1995  return ulDerivative;
1996 }
1997 
1998 
1999 scalar NURBS3DSurface::lengthDerivativeV
2000 (
2001  const scalar uConst,
2002  const scalar vStart,
2003  const scalar vEnd,
2004  const label nPts
2005 ) const
2006 {
2007  // Compute derivatives wrt v for the given v interval.
2008  vectorField dxdv(nPts, Zero);
2009  vectorField d2xdv2(nPts, Zero);
2010  scalarField localV(nPts, Zero);
2011  scalar vlDerivative(Zero);
2012 
2013  forAll(localV, vI)
2014  {
2015  scalar& vLocal(localV[vI]);
2016  vLocal = vStart + scalar(vI)/scalar(nPts-1)*(vEnd-vStart);
2017  dxdv[vI] = surfaceDerivativeV(uConst, vLocal);
2018  d2xdv2[vI] = surfaceDerivativeVV(uConst, vLocal);
2019  }
2020 
2021  // Integrate.
2022  for(label vI=0; vI<(nPts-1); vI++)
2023  {
2024  vlDerivative +=
2025  0.5
2026  * (
2027  (dxdv[vI+1]&d2xdv2[vI+1])/(mag(dxdv[vI+1])+SMALL)
2028  + (dxdv[vI ]&d2xdv2[vI ])/(mag(dxdv[vI ])+SMALL)
2029  )
2030  * (localV[vI+1]-localV[vI]);
2031  }
2032 
2033  return vlDerivative;
2034 }
2035 
2036 
2037 // * * * * * * * * * * * * * * * Access Functions * * * * * * * * * * * * * //
2038 
2039 const labelList& NURBS3DSurface::getBoundaryCPIDs()
2040 {
2041  if (!boundaryCPIDs_)
2042  {
2043  const label uNCPs(uBasis_.nCPs());
2044  const label vNCPs(vBasis_.nCPs());
2045  const label nBoundCPs(2*uNCPs+2*vNCPs-4);
2046  boundaryCPIDs_.reset(new labelList(nBoundCPs, -1));
2047  whichBoundaryCPID_.reset(new labelList(uNCPs*vNCPs, -1));
2048 
2049  // v-constant cps
2050  label bID(0);
2051  for(label vI=0; vI<vNCPs; vI+=vNCPs-1)
2052  {
2053  for(label uI=0; uI<uNCPs; uI++)
2054  {
2055  const label CPI(vI*uNCPs + uI);
2056  whichBoundaryCPID_()[CPI] = bID;
2057  boundaryCPIDs_()[bID++] = CPI;
2058  }
2059  }
2060  // u-constant cps
2061  for(label uI=0; uI<uNCPs; uI+=uNCPs-1)
2062  {
2063  // corner CPS already accounted for
2064  for(label vI=1; vI<vNCPs-1; vI++)
2065  {
2066  const label CPI(vI*uNCPs + uI);
2067  whichBoundaryCPID_()[CPI] = bID;
2068  boundaryCPIDs_()[bID++] = CPI;
2069  }
2070  }
2071  }
2072 
2073  return boundaryCPIDs_();
2074 }
2075 
2076 
2077 const labelList& NURBS3DSurface::getBoundaryCPIs()
2078 {
2079  return getBoundaryCPIDs();
2080 }
2081 
2082 
2083 const label& NURBS3DSurface::whichBoundaryCPI(const label& globalCPI)
2084 {
2085  if (!whichBoundaryCPID_)
2086  {
2087  getBoundaryCPIDs();
2088  }
2089 
2090  return whichBoundaryCPID_()[globalCPI];
2091 }
2092 
2093 
2094 // * * * * * * * * * * * * * * * Write Functions * * * * * * * * * * * * * //
2095 
2096 void NURBS3DSurface::write()
2097 {
2098  write(name_);
2099 }
2100 
2101 
2102 void NURBS3DSurface::write(const word fileName)
2103 {
2104  if (Pstream::master())
2105  {
2106  OFstream surfaceFile(fileName);
2107  OFstream surfaceFileCPs(fileName+"CPs");
2108  vectorField& surface(*this);
2109 
2110  forAll(*this, ptI)
2111  {
2112  surfaceFile
2113  << surface[ptI].x() << " "
2114  << surface[ptI].y() << " "
2115  << surface[ptI].z()
2116  << endl;
2117  }
2118 
2119  forAll(CPs_, CPI)
2120  {
2121  surfaceFileCPs
2122  << CPs_[CPI].x() << " "
2123  << CPs_[CPI].y() << " "
2124  << CPs_[CPI].z()
2125  << endl;
2126  }
2127 /*
2128  const label uDegree(uBasis_.degree());
2129  const label vDegree(vBasis_.degree());
2130  const label uNCPs(uBasis_.nCPs());
2131  const label vNCPs(vBasis_.nCPs());
2132 
2133  OFstream surfaceFileUBases(fileName+"UBases");
2134  OFstream surfaceFileVBases(fileName+"VBases");
2135 
2136  forAll(*this, ptI)
2137  {
2138  const scalar& u(u_[ptI]);
2139  const scalar& v(v_[ptI]);
2140  scalarField uBasesValues(uNCPs);
2141  scalarField vBasesValues(vNCPs);
2142 
2143  surfaceFileUBases << u << " ";
2144  surfaceFileVBases << v << " ";
2145 
2146  forAll(CPs_, CPI)
2147  {
2148  basesValues[CPI] = basis_.basisValue(CPI, degree, u);
2149  curveFileBases << basesValues[CPI] << " ";
2150  }
2151 
2152  curveFileBases << endl;
2153  }
2154 //*/
2155  }
2156 }
2157 
2158 
2159 void NURBS3DSurface::write(const fileName dirName, const fileName fileName)
2160 {
2161  if (Pstream::master())
2162  {
2163  OFstream surfaceFile(dirName/fileName);
2164  OFstream surfaceFileCPs(dirName/fileName+"CPs");
2165  vectorField& surface(*this);
2166 
2167  forAll(*this, ptI)
2168  {
2169  surfaceFile
2170  << surface[ptI].x() << " "
2171  << surface[ptI].y() << " "
2172  << surface[ptI].z()
2173  << endl;
2174  }
2175 
2176  forAll(CPs_, CPI)
2177  {
2178  surfaceFileCPs
2179  << CPs_[CPI].x() << " "
2180  << CPs_[CPI].y() << " "
2181  << CPs_[CPI].z()
2182  << endl;
2183  }
2184 /*
2185  const label uDegree(uBasis_.degree());
2186  const label vDegree(vBasis_.degree());
2187  const label uNCPs(uBasis_.nCPs());
2188  const label vNCPs(vBasis_.nCPs());
2189 
2190  OFstream surfaceFileUBases(dirName/fileName+"UBases");
2191  OFstream surfaceFileVBases(dirName/fileName+"VBases");
2192 
2193  forAll(*this, ptI)
2194  {
2195  const scalar& u(u_[ptI]);
2196  const scalar& v(v_[ptI]);
2197  scalarField uBasesValues(uNCPs);
2198  scalarField vBasesValues(vNCPs);
2199 
2200  surfaceFileUBases << u << " ";
2201  surfaceFileVBases << v << " ";
2202 
2203  forAll(CPs_, CPI)
2204  {
2205  basesValues[CPI] = basis_.basisValue(CPI, degree, u);
2206  curveFileBases << basesValues[CPI] << " ";
2207  }
2208 
2209  curveFileBases << endl;
2210  }
2211 //*/
2212  }
2213 }
2214 
2215 
2216 void NURBS3DSurface::writeWParses()
2217 {
2218  writeWParses(name_);
2219 }
2220 
2221 
2222 void NURBS3DSurface::writeWParses(const word fileName)
2223 {
2224  if (Pstream::master())
2225  {
2226  OFstream surfaceFile(fileName);
2227  OFstream surfaceFileCPs(fileName+"CPs");
2228  vectorField& surface(*this);
2229 
2230  forAll(*this, ptI)
2231  {
2232  surfaceFile
2233  << "("
2234  << surface[ptI].x() << " "
2235  << surface[ptI].y() << " "
2236  << surface[ptI].z() << ")"
2237  << endl;
2238  }
2239 
2240  forAll(CPs_, CPI)
2241  {
2242  surfaceFileCPs
2243  << "("
2244  << CPs_[CPI].x() << " "
2245  << CPs_[CPI].y() << " "
2246  << CPs_[CPI].z() << ")"
2247  << endl;
2248  }
2249 /*
2250  const label uDegree(uBasis_.degree());
2251  const label vDegree(vBasis_.degree());
2252  const label uNCPs(uBasis_.nCPs());
2253  const label vNCPs(vBasis_.nCPs());
2254 
2255  OFstream surfaceFileUBases(fileName+"UBases");
2256  OFstream surfaceFileVBases(fileName+"VBases");
2257 
2258  forAll(*this, ptI)
2259  {
2260  const scalar& u(u_[ptI]);
2261  const scalar& v(v_[ptI]);
2262  scalarField uBasesValues(uNCPs);
2263  scalarField vBasesValues(vNCPs);
2264 
2265  surfaceFileUBases << u << " ";
2266  surfaceFileVBases << v << " ";
2267 
2268  forAll(CPs_, CPI)
2269  {
2270  basesValues[CPI] = basis_.basisValue(CPI, degree, u);
2271  curveFileBases << basesValues[CPI] << " ";
2272  }
2273 
2274  curveFileBases << endl;
2275  }
2276 //*/
2277  }
2278 }
2279 
2280 
2281 void NURBS3DSurface::writeWParses
2282 (
2283  const fileName dirName,
2284  const fileName fileName
2285 )
2286 {
2287  if (Pstream::master())
2288  {
2289  OFstream surfaceFile(dirName/fileName);
2290  OFstream surfaceFileCPs(dirName/fileName+"CPs");
2291  vectorField& surface(*this);
2292 
2293  forAll(*this, ptI)
2294  {
2295  surfaceFile
2296  << "("
2297  << surface[ptI].x() << " "
2298  << surface[ptI].y() << " "
2299  << surface[ptI].z() << ")"
2300  << endl;
2301  }
2302 
2303  forAll(CPs_, CPI)
2304  {
2305  surfaceFileCPs
2306  << "("
2307  << CPs_[CPI].x() << " "
2308  << CPs_[CPI].y() << " "
2309  << CPs_[CPI].z() << ")"
2310  << endl;
2311  }
2312 /*
2313  const label uDegree(uBasis_.degree());
2314  const label vDegree(vBasis_.degree());
2315  const label uNCPs(uBasis_.nCPs());
2316  const label vNCPs(vBasis_.nCPs());
2317 
2318  OFstream surfaceFileUBases(dirName/fileName+"UBases");
2319  OFstream surfaceFileVBases(dirName/fileName+"VBases");
2320 
2321  forAll(*this, ptI)
2322  {
2323  const scalar& u(u_[ptI]);
2324  const scalar& v(v_[ptI]);
2325  scalarField uBasesValues(uNCPs);
2326  scalarField vBasesValues(vNCPs);
2327 
2328  surfaceFileUBases << u << " ";
2329  surfaceFileVBases << v << " ";
2330 
2331  forAll(CPs_, CPI)
2332  {
2333  basesValues[CPI] = basis_.basisValue(CPI, degree, u);
2334  curveFileBases << basesValues[CPI] << " ";
2335  }
2336 
2337  curveFileBases << endl;
2338  }
2339 */
2340  }
2341 }
2342 
2343 
2344 void NURBS3DSurface::writeVTK
2345 (
2346  const fileName vtkDirName,
2347  const fileName vtkFileName
2348 )
2349 {
2350  if (Pstream::master())
2351  {
2352  if (vtkFileName.ext() != word::null)
2353  {
2354  FatalErrorInFunction
2355  << "Do not supply a file extension."
2356  << exit(FatalError);
2357  }
2358 
2359  // Build the surface.
2360  buildSurface();
2361 
2362  // Write the surface as vtk.
2363  OFstream surfaceFile(vtkFileName);
2364  pointField& surfacePoints(*this);
2365  faceList surfaceFaces((nUPts_ - 1)*(nUPts_ - 1), face(4));
2366 
2367  for (label fuI = 0; fuI < (nUPts_ - 1); fuI++)
2368  {
2369  for (label fvI = 0; fvI < (nVPts_ - 1); fvI++)
2370  {
2371  const label fI(fuI*(nUPts_ - 1) + fvI);
2372  face& surfaceFace(surfaceFaces[fI]);
2373 
2374  surfaceFace[0] = (fuI)*nVPts_ + (fvI);
2375  surfaceFace[1] = (fuI + 1)*nVPts_ + (fvI);
2376  surfaceFace[2] = (fuI + 1)*nVPts_ + (fvI + 1);
2377  surfaceFace[3] = (fuI)*nVPts_ + (fvI + 1);
2378  }
2379  }
2380 
2381  surfaceWriters::vtkWriter writer;
2382 
2383  writer.open
2384  (
2385  surfacePoints,
2386  surfaceFaces,
2387  vtkDirName/vtkFileName,
2388  false
2389  );
2390 
2391  writer.close();
2392 
2393  // Write the control surface as vtk.
2394  fileName vtkCPFileName(vtkFileName+"CPs");
2395  pointField surfaceCPPoints(CPs_);
2396  const label uNCPs(uBasis_.nCPs());
2397  const label vNCPs(vBasis_.nCPs());
2398  faceList surfaceCPFaces((uNCPs-1)*(vNCPs-1), face(4));
2399 
2400  for (label fvCPI=0; fvCPI<(vNCPs-1); fvCPI++)
2401  {
2402  for (label fuCPI=0; fuCPI<(uNCPs-1); fuCPI++)
2403  {
2404  const label fCPI(fvCPI*(uNCPs-1) + fuCPI);
2405  face& surfaceCPFace(surfaceCPFaces[fCPI]);
2406 
2407  surfaceCPFace[0] = (fvCPI)*uNCPs + (fuCPI);
2408  surfaceCPFace[1] = (fvCPI + 1)*uNCPs + (fuCPI);
2409  surfaceCPFace[2] = (fvCPI + 1)*uNCPs + (fuCPI + 1);
2410  surfaceCPFace[3] = (fvCPI)*uNCPs + (fuCPI + 1);
2411  }
2412  }
2413 
2414  writer.open
2415  (
2416  surfaceCPPoints,
2417  surfaceCPFaces,
2418  vtkDirName/vtkCPFileName,
2419  false
2420  );
2421 
2422  writer.close();
2423  }
2424 }
2425 
2426 
2427 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2428 
2429 } // End namespace Foam
2430 
2431 // ************************************************************************* //
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:63
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:38
Foam::maxOp
Definition: ops.H:217
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:46
Foam::scalarList
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:61
Foam::NURBS3DSurface::surfacePoint
vector surfacePoint(const scalar &u, const scalar &v)
Definition: NURBS3DSurface.C:871
Foam::NURBSbasis::basisDerivativeU
scalar basisDerivativeU(const label iCP, const label degree, const scalar u) const
Definition: NURBSbasis.C:183
Foam::NURBS3DSurface::whichBoundaryCPI
const label & whichBoundaryCPI(const label &globalCPI)
Definition: NURBS3DSurface.C:2076
Foam::NURBS3DSurface::surfaceDerivativeU
vector surfaceDerivativeU(const scalar u, const scalar v) const
Definition: NURBS3DSurface.C:1612
Foam::NURBS3DSurface::surfaceDerivativeUU
vector surfaceDerivativeUU(const scalar u, const scalar v) const
Definition: NURBS3DSurface.C:1758
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:63
Foam::NURBS3DSurface::lengthDerivativeV
scalar lengthDerivativeV(const scalar uConst, const scalar vStart, const scalar vEnd, const label nPts) const
Definition: NURBS3DSurface.C:1993
Foam::fileName
A class for handling file names.
Definition: fileName.H:71
Foam::NURBS3DSurface::write
void write()
Definition: NURBS3DSurface.C:2089
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:57
Foam::Zero
static constexpr const zero Zero
Definition: zero.H:131
Foam::NURBS3DSurface::flipNrmOrientation
void flipNrmOrientation()
Definition: NURBS3DSurface.C:606
Foam::bound
volScalarField & bound(volScalarField &, const dimensionedScalar &lowerBound)
Definition: bound.C:28
Foam::NURBSbasis::basisValue
scalar basisValue(const label iCP, const label degree, const scalar u) const
Definition: NURBSbasis.C:133
Foam::NURBSbasis
NURBSbasis function. Used to construct NURBS curves, surfaces and volumes.
Definition: NURBSbasis.H:51
Foam::NURBS3DSurface::invertUV
void invertUV()
Definition: NURBS3DSurface.C:754
Foam::NURBS3DSurface::NURBS3DSurface
NURBS3DSurface(const List< vector > &CPs, const label nPointsU, const label nPointsV, const NURBSbasis &uBasis, const NURBSbasis &vBasis, const word name="NURBS3DSurface")
Definition: NURBS3DSurface.C:301
Foam::NURBS3DSurface::setName
void setName(const word &name)
Definition: NURBS3DSurface.C:631
Foam::NURBS3DSurface::writeVTK
void writeVTK(const fileName vtkDirName, const fileName vtkFileName)
Definition: NURBS3DSurface.C:2338
Foam::UPstream::master
static bool master(const label communicator=worldComm)
Definition: UPstream.H:453
Foam::NURBSbasis::degree
const label & degree() const
Definition: NURBSbasisI.H:33
Foam::endl
Ostream & endl(Ostream &os)
Definition: Ostream.H:381
Foam::NURBS3DSurface::writeWParses
void writeWParses()
Definition: NURBS3DSurface.C:2209
Foam::NURBS3DSurface::lengthDerivativeU
scalar lengthDerivativeU(const scalar vConst, const scalar uStart, const scalar uEnd, const label nPts) const
Definition: NURBS3DSurface.C:1955
Foam::NURBS3DSurface::checkRangeV
bool checkRangeV(const scalar v, const label CPI, const label vDegree) const
Definition: NURBS3DSurface.C:1399
Foam::sumOp
Definition: ops.H:207
Foam::NURBS3DSurface::surfaceDerivativeW
vector surfaceDerivativeW(const scalar u, const scalar v, const label CPI) const
Definition: NURBS3DSurface.C:1912
forAll
#define forAll(list, i)
Definition: stdFoam.H:349
OFstream.H
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:48
Foam::NURBS3DSurface::checkRangeUV
bool checkRangeUV(const scalar v, const scalar u, const label CPI, const label uDegree, const label vDegree) const
Definition: NURBS3DSurface.C:1424
Foam::NURBSbasis::basisDerivativeUU
scalar basisDerivativeUU(const label iCP, const label degree, const scalar u) const
Definition: NURBSbasis.C:224
Foam::NURBS3DSurface::checkRangeU
bool checkRangeU(const scalar u, const label CPI, const label uDegree) const
Definition: NURBS3DSurface.C:1374
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:45
R
#define R(A, B, C, D, E, F, K, M)
Foam::surfaceWriters::vtkWriter
A surfaceWriter for VTK legacy (.vtk) or XML (.vtp) format.
Definition: vtkSurfaceWriter.H:126
Foam::constant::physicoChemical::b
const dimensionedScalar b
Definition: createFields.H:27
Foam::Field
Generic templated field type.
Definition: Field.H:59
Foam::Info
messageStream Info
Foam::List::setSize
void setSize(const label n)
Definition: List.H:218
Foam::NURBS3DSurface::invertV
void invertV()
Definition: NURBS3DSurface.C:725
delta
scalar delta
Definition: LISASMDCalcMethod2.H:8
Foam::NURBS3DSurface::surfaceDerivativeV
vector surfaceDerivativeV(const scalar u, const scalar v) const
Definition: NURBS3DSurface.C:1654
Foam::NURBS3DSurface::findClosestSurfacePoint
scalarList findClosestSurfacePoint(const vector &targetPoint, const label maxIter=100, const scalar tolerance=1.e-6)
Definition: NURBS3DSurface.C:918
Foam::NURBS3DSurface::setNrmOrientation
void setNrmOrientation(const vector &givenNrm, const scalar u, const scalar v)
Definition: NURBS3DSurface.C:578
field
rDeltaTY field()
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Definition: hashSets.C:40
Foam::FatalError
error FatalError
Foam::NURBS3DSurface::lengthV
scalar lengthV(const label uIConst, const label vIStart, const label vIEnd) const
Definition: NURBS3DSurface.C:1533
Foam::writer
Base class for graphics format writing. Entry points are.
Definition: writer.H:77
Foam::NURBSbasis::nCPs
const label & nCPs() const
Definition: NURBSbasisI.H:39
Foam::fileName::ext
word ext() const
Definition: fileNameI.H:211
Foam::NURBS3DSurface::setWeights
void setWeights(const scalarList &weights)
Definition: NURBS3DSurface.C:625
Foam
Definition: atmBoundaryLayer.C:26
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:47
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:125
Foam::NURBS3DSurface::getBoundaryCPIs
const labelList & getBoundaryCPIs()
Definition: NURBS3DSurface.C:2070
Foam::NURBS3DSurface::invertU
void invertU()
Definition: NURBS3DSurface.C:696
Foam::OFstream
Output to file stream, using an OSstream.
Definition: OFstream.H:49
U
U
Definition: pEqn.H:72
Foam::NURBS3DSurface::surfaceDerivativeUV
vector surfaceDerivativeUV(const scalar u, const scalar v) const
Definition: NURBS3DSurface.C:1696
FatalErrorInFunction
#define FatalErrorInFunction
Definition: error.H:465
Foam::faceList
List< face > faceList
A List of faces.
Definition: faceListFwd.H:41
Foam::Vector< scalar >
Foam::NURBS3DSurface::nrm
const vector nrm(scalar u, scalar v)
Definition: NURBS3DSurface.C:1258
Foam::NURBS3DSurface::lengthU
scalar lengthU(const label vIConst, const label uIStart, const label uIEnd) const
Definition: NURBS3DSurface.C:1456
vtkSurfaceWriter.H
Foam::List< vector >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::constant::electromagnetic::e
const dimensionedScalar e
Definition: createFields.H:11
Foam::tmp::New
static tmp< T > New(Args &&... args)
Foam::word::null
static const word null
Definition: word.H:78
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::constant::universal::c
const dimensionedScalar c
Foam::NURBS3DSurface::surfaceDerivativeVV
vector surfaceDerivativeVV(const scalar u, const scalar v) const
Definition: NURBS3DSurface.C:1814
Foam::name
word name(const expressions::valueTypeCode typeCode)
Definition: exprTraits.C:52
NURBS3DSurface.H
Foam::VectorSpace< Vector< scalar >, scalar, 3 >::zero
static const Form zero
Definition: VectorSpace.H:111
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:68
Foam::NURBS3DSurface::makeEquidistant
void makeEquidistant(const label lenAcc=25, const label maxIter=10, const label spacingCorrInterval=-1, const scalar tolerance=1.e-5)
Definition: NURBS3DSurface.C:785
Foam::point
vector point
Point is a vector.
Definition: point.H:37
writer
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
Foam::NURBS3DSurface::getBoundaryCPIDs
const labelList & getBoundaryCPIDs()
Definition: NURBS3DSurface.C:2032
Foam::NURBS3DSurface::genEquidistant
List< scalarList > genEquidistant(const label nUPts=100, const label nVPts=100, const label lenAcc=25, const label maxIter=10, const label spacingCorrInterval=-1, const scalar tolerance=1.e-5)
Definition: NURBS3DSurface.C:1278
WarningInFunction
#define WarningInFunction
Definition: messageStream.H:353
Foam::NURBSbasis::checkRange
bool checkRange(const scalar u, const label CPI, const label degree) const
Definition: NURBSbasis.C:266
Foam::NURBS3DSurface::buildSurface
void buildSurface()
Definition: NURBS3DSurface.C:637
Foam::NURBS3DSurface::surfaceDerivativeCP
scalar surfaceDerivativeCP(const scalar u, const scalar v, const label CPI) const
Definition: NURBS3DSurface.C:1870
Foam::NURBS3DSurface::setCPs
void setCPs(const List< vector > &CPs)
Definition: NURBS3DSurface.C:619

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17