MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
nurbs.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "mesh_headers.hpp"
13 #include "../fem/fem.hpp"
14 #include "../general/text.hpp"
15 
16 #include <fstream>
17 #include <algorithm>
18 #if defined(_MSC_VER) && (_MSC_VER < 1800)
19 #include <float.h>
20 #define copysign _copysign
21 #endif
22 
23 namespace mfem
24 {
25 
26 using namespace std;
27 
28 const int KnotVector::MaxOrder = 10;
29 
30 KnotVector::KnotVector(std::istream &input)
31 {
32  input >> Order >> NumOfControlPoints;
33  knot.Load(input, NumOfControlPoints + Order + 1);
34  GetElements();
35 }
36 
37 KnotVector::KnotVector(int Order_, int NCP)
38 {
39  Order = Order_;
40  NumOfControlPoints = NCP;
41  knot.SetSize(NumOfControlPoints + Order + 1);
42 
43  knot = -1.;
44 }
45 
47 {
48  Order = kv.Order;
49  NumOfControlPoints = kv.NumOfControlPoints;
50  NumOfElements = kv.NumOfElements;
51  knot = kv.knot;
52  // alternatively, re-compute NumOfElements
53  // GetElements();
54 
55  return *this;
56 }
57 
59 {
60  if (t < 0)
61  {
62  mfem_error("KnotVector::DegreeElevate :\n"
63  " Parent KnotVector order higher than child");
64  }
65 
66  const int nOrder = Order + t;
67  KnotVector *newkv = new KnotVector(nOrder, GetNCP() + t);
68 
69  for (int i = 0; i <= nOrder; i++)
70  {
71  (*newkv)[i] = knot(0);
72  }
73  for (int i = nOrder + 1; i < newkv->GetNCP(); i++)
74  {
75  (*newkv)[i] = knot(i - t);
76  }
77  for (int i = 0; i <= nOrder; i++)
78  {
79  (*newkv)[newkv->GetNCP() + i] = knot(knot.Size()-1);
80  }
81 
82  newkv->GetElements();
83 
84  return newkv;
85 }
86 
88 {
89  newknots.SetSize(NumOfElements);
90  int j = 0;
91  for (int i = 0; i < knot.Size()-1; i++)
92  {
93  if (knot(i) != knot(i+1))
94  {
95  newknots(j) = 0.5*(knot(i) + knot(i+1));
96  j++;
97  }
98  }
99 }
100 
102 {
103  NumOfElements = 0;
104  for (int i = Order; i < NumOfControlPoints; i++)
105  {
106  if (knot(i) != knot(i+1))
107  {
108  NumOfElements++;
109  }
110  }
111 }
112 
114 {
115  double apb = knot(0) + knot(knot.Size()-1);
116 
117  int ns = (NumOfControlPoints - Order)/2;
118  for (int i = 1; i <= ns; i++)
119  {
120  double tmp = apb - knot(Order + i);
121  knot(Order + i) = apb - knot(NumOfControlPoints - i);
122  knot(NumOfControlPoints - i) = tmp;
123  }
124 }
125 
126 void KnotVector::Print(std::ostream &os) const
127 {
128  os << Order << ' ' << NumOfControlPoints << ' ';
129  knot.Print(os, knot.Size());
130 }
131 
132 
133 void KnotVector::PrintFunctions(std::ostream &os, int samples) const
134 {
135  Vector shape(Order+1);
136 
137  double x, dx = 1.0/double (samples - 1);
138 
139  for (int i = 0; i <GetNE() ; i++)
140  {
141  for (int j = 0; j <samples; j++)
142  {
143  x =j*dx;
144  os<< x + i;
145 
146  CalcShape ( shape, i, x);
147  for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
148 
149  CalcDShape ( shape, i, x);
150  for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
151 
152  CalcD2Shape ( shape, i, x);
153  for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
154  os<<endl;
155  }
156  }
157 }
158 
159 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
160 void KnotVector::CalcShape(Vector &shape, int i, double xi) const
161 {
162  MFEM_ASSERT(Order <= MaxOrder, "Order > MaxOrder!");
163 
164  int p = Order;
165  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
166  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), saved, tmp;
167  double left[MaxOrder+1], right[MaxOrder+1];
168 
169  shape(0) = 1.;
170  for (int j = 1; j <= p; ++j)
171  {
172  left[j] = u - knot(ip+1-j);
173  right[j] = knot(ip+j) - u;
174  saved = 0.;
175  for (int r = 0; r < j; ++r)
176  {
177  tmp = shape(r)/(right[r+1] + left[j-r]);
178  shape(r) = saved + right[r+1]*tmp;
179  saved = left[j-r]*tmp;
180  }
181  shape(j) = saved;
182  }
183 }
184 
185 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
186 void KnotVector::CalcDShape(Vector &grad, int i, double xi) const
187 {
188  int p = Order, rk, pk;
189  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
190  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), temp, saved, d;
191  double ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1], right[MaxOrder+1];
192 
193 #ifdef MFEM_DEBUG
194  if (p > MaxOrder)
195  {
196  mfem_error("KnotVector::CalcDShape : Order > MaxOrder!");
197  }
198 #endif
199 
200  ndu[0][0] = 1.0;
201  for (int j = 1; j <= p; j++)
202  {
203  left[j] = u - knot(ip-j+1);
204  right[j] = knot(ip+j) - u;
205  saved = 0.0;
206  for (int r = 0; r < j; r++)
207  {
208  ndu[j][r] = right[r+1] + left[j-r];
209  temp = ndu[r][j-1]/ndu[j][r];
210  ndu[r][j] = saved + right[r+1]*temp;
211  saved = left[j-r]*temp;
212  }
213  ndu[j][j] = saved;
214  }
215 
216  for (int r = 0; r <= p; ++r)
217  {
218  d = 0.0;
219  rk = r-1;
220  pk = p-1;
221  if (r >= 1)
222  {
223  d = ndu[rk][pk]/ndu[p][rk];
224  }
225  if (r <= pk)
226  {
227  d -= ndu[r][pk]/ndu[p][r];
228  }
229  grad(r) = d;
230  }
231 
232  if (i >= 0)
233  {
234  grad *= p*(knot(ip+1) - knot(ip));
235  }
236  else
237  {
238  grad *= p*(knot(ip) - knot(ip+1));
239  }
240 }
241 
242 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
243 void KnotVector::CalcDnShape(Vector &gradn, int n, int i, double xi) const
244 {
245  int p = Order, rk, pk, j1, j2,r,j,k;
246  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
247  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip);
248  double temp, saved, d;
249  double a[2][MaxOrder+1],ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1],
250  right[MaxOrder+1];
251 
252 #ifdef MFEM_DEBUG
253  if (p > MaxOrder)
254  {
255  mfem_error("KnotVector::CalcDnShape : Order > MaxOrder!");
256  }
257 #endif
258 
259  ndu[0][0] = 1.0;
260  for (j = 1; j <= p; j++)
261  {
262  left[j] = u - knot(ip-j+1);
263  right[j] = knot(ip+j)- u;
264 
265  saved = 0.0;
266  for (r = 0; r < j; r++)
267  {
268  ndu[j][r] = right[r+1] + left[j-r];
269  temp = ndu[r][j-1]/ndu[j][r];
270  ndu[r][j] = saved + right[r+1]*temp;
271  saved = left[j-r]*temp;
272  }
273  ndu[j][j] = saved;
274  }
275 
276  for (r = 0; r <= p; r++)
277  {
278  int s1 = 0;
279  int s2 = 1;
280  a[0][0] = 1.0;
281  for (k = 1; k <= n; k++)
282  {
283  d = 0.0;
284  rk = r-k;
285  pk = p-k;
286  if (r >= k)
287  {
288  a[s2][0] = a[s1][0]/ndu[pk+1][rk];
289  d = a[s2][0]*ndu[rk][pk];
290  }
291 
292  if (rk >= -1)
293  {
294  j1 = 1;
295  }
296  else
297  {
298  j1 = -rk;
299  }
300 
301  if (r-1<= pk)
302  {
303  j2 = k-1;
304  }
305  else
306  {
307  j2 = p-r;
308  }
309 
310  for (j = j1; j <= j2; j++)
311  {
312  a[s2][j] = (a[s1][j] - a[s1][j-1])/ndu[pk+1][rk+j];
313  d += a[s2][j]*ndu[rk+j][pk];
314  }
315 
316  if (r <= pk)
317  {
318  a[s2][k] = - a[s1][k-1]/ndu[pk+1][r];
319  d += a[s2][j]*ndu[rk+j][pk];
320  }
321  gradn[r] = d;
322  j = s1;
323  s1 = s2;
324  s2 = j;
325  }
326  }
327 
328  if (i >= 0)
329  {
330  u = (knot(ip+1) - knot(ip));
331  }
332  else
333  {
334  u = (knot(ip) - knot(ip+1));
335  }
336 
337  temp = p*u;
338  for (k = 1; k <= n-1; k++) { temp *= (p-k)*u; }
339 
340  for (j = 0; j <= p; j++) { gradn[j] *= temp; }
341 
342 }
343 
344 
345 int KnotVector::findKnotSpan(double u) const
346 {
347  int low, mid, high;
348 
349  if (u == knot(NumOfControlPoints+Order))
350  {
351  mid = NumOfControlPoints;
352  }
353  else
354  {
355  low = Order;
356  high = NumOfControlPoints + 1;
357  mid = (low + high)/2;
358  while ( (u < knot(mid-1)) || (u > knot(mid)) )
359  {
360  if (u < knot(mid-1))
361  {
362  high = mid;
363  }
364  else
365  {
366  low = mid;
367  }
368  mid = (low + high)/2;
369  }
370  }
371  return mid;
372 }
373 
374 void KnotVector::Difference(const KnotVector &kv, Vector &diff) const
375 {
376  if (Order != kv.GetOrder())
377  {
378  mfem_error("KnotVector::Difference :\n"
379  " Can not compare knot vectors with different orders!");
380  }
381 
382  int s = kv.Size() - Size();
383  if (s < 0)
384  {
385  kv.Difference(*this, diff);
386  return;
387  }
388 
389  diff.SetSize(s);
390 
391  s = 0;
392  int i = 0;
393  for (int j = 0; j < kv.Size(); j++)
394  {
395  if (knot(i) == kv[j])
396  {
397  i++;
398  }
399  else
400  {
401  diff(s) = kv[j];
402  s++;
403  }
404  }
405 }
406 
407 void NURBSPatch::init(int dim_)
408 {
409  Dim = dim_;
410  sd = nd = -1;
411 
412  if (kv.Size() == 2)
413  {
414  ni = kv[0]->GetNCP();
415  nj = kv[1]->GetNCP();
416  nk = -1;
417 
418  data = new double[ni*nj*Dim];
419 
420 #ifdef MFEM_DEBUG
421  for (int i = 0; i < ni*nj*Dim; i++)
422  {
423  data[i] = -999.99;
424  }
425 #endif
426  }
427  else if (kv.Size() == 3)
428  {
429  ni = kv[0]->GetNCP();
430  nj = kv[1]->GetNCP();
431  nk = kv[2]->GetNCP();
432 
433  data = new double[ni*nj*nk*Dim];
434 
435 #ifdef MFEM_DEBUG
436  for (int i = 0; i < ni*nj*nk*Dim; i++)
437  {
438  data[i] = -999.99;
439  }
440 #endif
441  }
442  else
443  {
444  mfem_error("NURBSPatch::init : Wrong dimension of knotvectors!");
445  }
446 }
447 
448 NURBSPatch::NURBSPatch(const NURBSPatch &orig)
449  : ni(orig.ni), nj(orig.nj), nk(orig.nk), Dim(orig.Dim),
450  data(NULL), kv(orig.kv.Size()), sd(orig.sd), nd(orig.nd)
451 {
452  // Allocate and copy data:
453  const int data_size = Dim*ni*nj*((kv.Size() == 2) ? 1 : nk);
454  data = new double[data_size];
455  std::memcpy(data, orig.data, data_size*sizeof(double));
456 
457  // Copy the knot vectors:
458  for (int i = 0; i < kv.Size(); i++)
459  {
460  kv[i] = new KnotVector(*orig.kv[i]);
461  }
462 }
463 
464 NURBSPatch::NURBSPatch(std::istream &input)
465 {
466  int pdim, dim, size = 1;
467  string ident;
468 
469  input >> ws >> ident >> pdim; // knotvectors
470  kv.SetSize(pdim);
471  for (int i = 0; i < pdim; i++)
472  {
473  kv[i] = new KnotVector(input);
474  size *= kv[i]->GetNCP();
475  }
476 
477  input >> ws >> ident >> dim; // dimension
478  init(dim + 1);
479 
480  input >> ws >> ident; // controlpoints (homogeneous coordinates)
481  if (ident == "controlpoints" || ident == "controlpoints_homogeneous")
482  {
483  for (int j = 0, i = 0; i < size; i++)
484  for (int d = 0; d <= dim; d++, j++)
485  {
486  input >> data[j];
487  }
488  }
489  else // "controlpoints_cartesian" (Cartesian coordinates with weight)
490  {
491  for (int j = 0, i = 0; i < size; i++)
492  {
493  for (int d = 0; d <= dim; d++)
494  {
495  input >> data[j+d];
496  }
497  for (int d = 0; d < dim; d++)
498  {
499  data[j+d] *= data[j+dim];
500  }
501  j += (dim+1);
502  }
503  }
504 }
505 
506 NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim_)
507 {
508  kv.SetSize(2);
509  kv[0] = new KnotVector(*kv0);
510  kv[1] = new KnotVector(*kv1);
511  init(dim_);
512 }
513 
514 NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
515  const KnotVector *kv2, int dim_)
516 {
517  kv.SetSize(3);
518  kv[0] = new KnotVector(*kv0);
519  kv[1] = new KnotVector(*kv1);
520  kv[2] = new KnotVector(*kv2);
521  init(dim_);
522 }
523 
524 NURBSPatch::NURBSPatch(Array<const KnotVector *> &kv_, int dim_)
525 {
526  kv.SetSize(kv_.Size());
527  for (int i = 0; i < kv.Size(); i++)
528  {
529  kv[i] = new KnotVector(*kv_[i]);
530  }
531  init(dim_);
532 }
533 
534 NURBSPatch::NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP)
535 {
536  kv.SetSize(parent->kv.Size());
537  for (int i = 0; i < kv.Size(); i++)
538  if (i != dir)
539  {
540  kv[i] = new KnotVector(*parent->kv[i]);
541  }
542  else
543  {
544  kv[i] = new KnotVector(Order, NCP);
545  }
546  init(parent->Dim);
547 }
548 
549 void NURBSPatch::swap(NURBSPatch *np)
550 {
551  if (data != NULL)
552  {
553  delete [] data;
554  }
555 
556  for (int i = 0; i < kv.Size(); i++)
557  {
558  if (kv[i]) { delete kv[i]; }
559  }
560 
561  data = np->data;
562  np->kv.Copy(kv);
563 
564  ni = np->ni;
565  nj = np->nj;
566  nk = np->nk;
567  Dim = np->Dim;
568 
569  np->data = NULL;
570  np->kv.SetSize(0);
571 
572  delete np;
573 }
574 
575 NURBSPatch::~NURBSPatch()
576 {
577  if (data != NULL)
578  {
579  delete [] data;
580  }
581 
582  for (int i = 0; i < kv.Size(); i++)
583  {
584  if (kv[i]) { delete kv[i]; }
585  }
586 }
587 
588 void NURBSPatch::Print(std::ostream &os) const
589 {
590  int size = 1;
591 
592  os << "knotvectors\n" << kv.Size() << '\n';
593  for (int i = 0; i < kv.Size(); i++)
594  {
595  kv[i]->Print(os);
596  size *= kv[i]->GetNCP();
597  }
598 
599  os << "\ndimension\n" << Dim - 1
600  << "\n\ncontrolpoints\n";
601  for (int j = 0, i = 0; i < size; i++)
602  {
603  os << data[j++];
604  for (int d = 1; d < Dim; d++)
605  {
606  os << ' ' << data[j++];
607  }
608  os << '\n';
609  }
610 }
611 
612 int NURBSPatch::SetLoopDirection(int dir)
613 {
614  if (nk == -1)
615  {
616  if (dir == 0)
617  {
618  sd = Dim;
619  nd = ni;
620 
621  return nj*Dim;
622  }
623  else if (dir == 1)
624  {
625  sd = ni*Dim;
626  nd = nj;
627 
628  return ni*Dim;
629  }
630  else
631  {
632  mfem::err << "NURBSPatch::SetLoopDirection :\n"
633  " Direction error in 2D patch, dir = " << dir << '\n';
634  mfem_error();
635  }
636  }
637  else
638  {
639  if (dir == 0)
640  {
641  sd = Dim;
642  nd = ni;
643 
644  return nj*nk*Dim;
645  }
646  else if (dir == 1)
647  {
648  sd = ni*Dim;
649  nd = nj;
650 
651  return ni*nk*Dim;
652  }
653  else if (dir == 2)
654  {
655  sd = ni*nj*Dim;
656  nd = nk;
657 
658  return ni*nj*Dim;
659  }
660  else
661  {
662  mfem::err << "NURBSPatch::SetLoopDirection :\n"
663  " Direction error in 3D patch, dir = " << dir << '\n';
664  mfem_error();
665  }
666  }
667 
668  return -1;
669 }
670 
671 void NURBSPatch::UniformRefinement()
672 {
673  Vector newknots;
674  for (int dir = 0; dir < kv.Size(); dir++)
675  {
676  kv[dir]->UniformRefinement(newknots);
677  KnotInsert(dir, newknots);
678  }
679 }
680 
681 void NURBSPatch::KnotInsert(Array<KnotVector *> &newkv)
682 {
683  for (int dir = 0; dir < kv.Size(); dir++)
684  {
685  KnotInsert(dir, *newkv[dir]);
686  }
687 }
688 
689 void NURBSPatch::KnotInsert(int dir, const KnotVector &newkv)
690 {
691  if (dir >= kv.Size() || dir < 0)
692  {
693  mfem_error("NURBSPatch::KnotInsert : Incorrect direction!");
694  }
695 
696  int t = newkv.GetOrder() - kv[dir]->GetOrder();
697 
698  if (t > 0)
699  {
700  DegreeElevate(dir, t);
701  }
702  else if (t < 0)
703  {
704  mfem_error("NURBSPatch::KnotInsert : Incorrect order!");
705  }
706 
707  Vector diff;
708  GetKV(dir)->Difference(newkv, diff);
709  if (diff.Size() > 0)
710  {
711  KnotInsert(dir, diff);
712  }
713 }
714 
715 void NURBSPatch::KnotInsert(Array<Vector *> &newkv)
716 {
717  for (int dir = 0; dir < kv.Size(); dir++)
718  {
719  KnotInsert(dir, *newkv[dir]);
720  }
721 }
722 
723 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
724 void NURBSPatch::KnotInsert(int dir, const Vector &knot)
725 {
726  if (knot.Size() == 0 ) { return; }
727 
728  if (dir >= kv.Size() || dir < 0)
729  {
730  mfem_error("NURBSPatch::KnotInsert : Incorrect direction!");
731  }
732 
733  NURBSPatch &oldp = *this;
734  KnotVector &oldkv = *kv[dir];
735 
736  NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder(),
737  oldkv.GetNCP() + knot.Size());
738  NURBSPatch &newp = *newpatch;
739  KnotVector &newkv = *newp.GetKV(dir);
740 
741  int size = oldp.SetLoopDirection(dir);
742  if (size != newp.SetLoopDirection(dir))
743  {
744  mfem_error("NURBSPatch::KnotInsert : Size mismatch!");
745  }
746 
747  int rr = knot.Size() - 1;
748  int a = oldkv.findKnotSpan(knot(0)) - 1;
749  int b = oldkv.findKnotSpan(knot(rr)) - 1;
750  int pl = oldkv.GetOrder();
751  int ml = oldkv.GetNCP();
752 
753  for (int j = 0; j <= a; j++)
754  {
755  newkv[j] = oldkv[j];
756  }
757  for (int j = b+pl; j <= ml+pl; j++)
758  {
759  newkv[j+rr+1] = oldkv[j];
760  }
761  for (int k = 0; k <= (a-pl); k++)
762  {
763  for (int ll = 0; ll < size; ll++)
764  {
765  newp(k,ll) = oldp(k,ll);
766  }
767  }
768  for (int k = (b-1); k < ml; k++)
769  {
770  for (int ll = 0; ll < size; ll++)
771  {
772  newp(k+rr+1,ll) = oldp(k,ll);
773  }
774  }
775 
776  int i = b+pl-1;
777  int k = b+pl+rr;
778 
779  for (int j = rr; j >= 0; j--)
780  {
781  while ( (knot(j) <= oldkv[i]) && (i > a) )
782  {
783  newkv[k] = oldkv[i];
784  for (int ll = 0; ll < size; ll++)
785  {
786  newp(k-pl-1,ll) = oldp(i-pl-1,ll);
787  }
788 
789  k--;
790  i--;
791  }
792 
793  for (int ll = 0; ll < size; ll++)
794  {
795  newp(k-pl-1,ll) = newp(k-pl,ll);
796  }
797 
798  for (int l = 1; l <= pl; l++)
799  {
800  int ind = k-pl+l;
801  double alfa = newkv[k+l] - knot(j);
802  if (fabs(alfa) == 0.0)
803  {
804  for (int ll = 0; ll < size; ll++)
805  {
806  newp(ind-1,ll) = newp(ind,ll);
807  }
808  }
809  else
810  {
811  alfa = alfa/(newkv[k+l] - oldkv[i-pl+l]);
812  for (int ll = 0; ll < size; ll++)
813  {
814  newp(ind-1,ll) = alfa*newp(ind-1,ll) + (1.0-alfa)*newp(ind,ll);
815  }
816  }
817  }
818 
819  newkv[k] = knot(j);
820  k--;
821  }
822 
823  newkv.GetElements();
824 
825  swap(newpatch);
826 }
827 
828 void NURBSPatch::DegreeElevate(int t)
829 {
830  for (int dir = 0; dir < kv.Size(); dir++)
831  {
832  DegreeElevate(dir, t);
833  }
834 }
835 
836 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
837 void NURBSPatch::DegreeElevate(int dir, int t)
838 {
839  if (dir >= kv.Size() || dir < 0)
840  {
841  mfem_error("NURBSPatch::DegreeElevate : Incorrect direction!");
842  }
843 
844  int i, j, k, kj, mpi, mul, mh, kind, cind, first, last;
845  int r, a, b, oldr, save, s, tr, lbz, rbz, l;
846  double inv, ua, ub, numer, alf, den, bet, gam;
847 
848  NURBSPatch &oldp = *this;
849  KnotVector &oldkv = *kv[dir];
850 
851  NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder() + t,
852  oldkv.GetNCP() + oldkv.GetNE()*t);
853  NURBSPatch &newp = *newpatch;
854  KnotVector &newkv = *newp.GetKV(dir);
855 
856  int size = oldp.SetLoopDirection(dir);
857  if (size != newp.SetLoopDirection(dir))
858  {
859  mfem_error("NURBSPatch::DegreeElevate : Size mismatch!");
860  }
861 
862  int p = oldkv.GetOrder();
863  int n = oldkv.GetNCP()-1;
864 
865  DenseMatrix bezalfs (p+t+1, p+1);
866  DenseMatrix bpts (p+1, size);
867  DenseMatrix ebpts (p+t+1, size);
868  DenseMatrix nextbpts(p-1, size);
869  Vector alphas (p-1);
870 
871  int m = n + p + 1;
872  int ph = p + t;
873  int ph2 = ph/2;
874 
875  {
876  Array2D<int> binom(ph+1, ph+1);
877  for (i = 0; i <= ph; i++)
878  {
879  binom(i,0) = binom(i,i) = 1;
880  for (j = 1; j < i; j++)
881  {
882  binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
883  }
884  }
885 
886  bezalfs(0,0) = 1.0;
887  bezalfs(ph,p) = 1.0;
888 
889  for (i = 1; i <= ph2; i++)
890  {
891  inv = 1.0/binom(ph,i);
892  mpi = min(p,i);
893  for (j = max(0,i-t); j <= mpi; j++)
894  {
895  bezalfs(i,j) = inv*binom(p,j)*binom(t,i-j);
896  }
897  }
898  }
899 
900  for (i = ph2+1; i < ph; i++)
901  {
902  mpi = min(p,i);
903  for (j = max(0,i-t); j <= mpi; j++)
904  {
905  bezalfs(i,j) = bezalfs(ph-i,p-j);
906  }
907  }
908 
909  mh = ph;
910  kind = ph + 1;
911  r = -1;
912  a = p;
913  b = p + 1;
914  cind = 1;
915  ua = oldkv[0];
916  for (l = 0; l < size; l++)
917  {
918  newp(0,l) = oldp(0,l);
919  }
920  for (i = 0; i <= ph; i++)
921  {
922  newkv[i] = ua;
923  }
924 
925  for (i = 0; i <= p; i++)
926  {
927  for (l = 0; l < size; l++)
928  {
929  bpts(i,l) = oldp(i,l);
930  }
931  }
932 
933  while (b < m)
934  {
935  i = b;
936  while (b < m && oldkv[b] == oldkv[b+1]) { b++; }
937 
938  mul = b-i+1;
939 
940  mh = mh + mul + t;
941  ub = oldkv[b];
942  oldr = r;
943  r = p-mul;
944  if (oldr > 0) { lbz = (oldr+2)/2; }
945  else { lbz = 1; }
946 
947  if (r > 0) { rbz = ph-(r+1)/2; }
948  else { rbz = ph; }
949 
950  if (r > 0)
951  {
952  numer = ub - ua;
953  for (k = p ; k > mul; k--)
954  {
955  alphas[k-mul-1] = numer/(oldkv[a+k]-ua);
956  }
957 
958  for (j = 1; j <= r; j++)
959  {
960  save = r-j;
961  s = mul+j;
962  for (k = p; k >= s; k--)
963  {
964  for (l = 0; l < size; l++)
965  bpts(k,l) = (alphas[k-s]*bpts(k,l) +
966  (1.0-alphas[k-s])*bpts(k-1,l));
967  }
968  for (l = 0; l < size; l++)
969  {
970  nextbpts(save,l) = bpts(p,l);
971  }
972  }
973  }
974 
975  for (i = lbz; i <= ph; i++)
976  {
977  for (l = 0; l < size; l++)
978  {
979  ebpts(i,l) = 0.0;
980  }
981  mpi = min(p,i);
982  for (j = max(0,i-t); j <= mpi; j++)
983  {
984  for (l = 0; l < size; l++)
985  {
986  ebpts(i,l) += bezalfs(i,j)*bpts(j,l);
987  }
988  }
989  }
990 
991  if (oldr > 1)
992  {
993  first = kind-2;
994  last = kind;
995  den = ub-ua;
996  bet = (ub-newkv[kind-1])/den;
997 
998  for (tr = 1; tr < oldr; tr++)
999  {
1000  i = first;
1001  j = last;
1002  kj = j-kind+1;
1003  while (j-i > tr)
1004  {
1005  if (i < cind)
1006  {
1007  alf = (ub-newkv[i])/(ua-newkv[i]);
1008  for (l = 0; l < size; l++)
1009  {
1010  newp(i,l) = alf*newp(i,l)-(1.0-alf)*newp(i-1,l);
1011  }
1012  }
1013  if (j >= lbz)
1014  {
1015  if ((j-tr) <= (kind-ph+oldr))
1016  {
1017  gam = (ub-newkv[j-tr])/den;
1018  for (l = 0; l < size; l++)
1019  {
1020  ebpts(kj,l) = gam*ebpts(kj,l) + (1.0-gam)*ebpts(kj+1,l);
1021  }
1022  }
1023  else
1024  {
1025  for (l = 0; l < size; l++)
1026  {
1027  ebpts(kj,l) = bet*ebpts(kj,l) + (1.0-bet)*ebpts(kj+1,l);
1028  }
1029  }
1030  }
1031  i = i+1;
1032  j = j-1;
1033  kj = kj-1;
1034  }
1035  first--;
1036  last++;
1037  }
1038  }
1039 
1040  if (a != p)
1041  {
1042  for (i = 0; i < (ph-oldr); i++)
1043  {
1044  newkv[kind] = ua;
1045  kind = kind+1;
1046  }
1047  }
1048  for (j = lbz; j <= rbz; j++)
1049  {
1050  for (l = 0; l < size; l++)
1051  {
1052  newp(cind,l) = ebpts(j,l);
1053  }
1054  cind = cind +1;
1055  }
1056 
1057  if (b < m)
1058  {
1059  for (j = 0; j <r; j++)
1060  for (l = 0; l < size; l++)
1061  {
1062  bpts(j,l) = nextbpts(j,l);
1063  }
1064 
1065  for (j = r; j <= p; j++)
1066  for (l = 0; l < size; l++)
1067  {
1068  bpts(j,l) = oldp(b-p+j,l);
1069  }
1070 
1071  a = b;
1072  b = b+1;
1073  ua = ub;
1074  }
1075  else
1076  {
1077  for (i = 0; i <= ph; i++)
1078  {
1079  newkv[kind+i] = ub;
1080  }
1081  }
1082  }
1083  newkv.GetElements();
1084 
1085  swap(newpatch);
1086 }
1087 
1088 void NURBSPatch::FlipDirection(int dir)
1089 {
1090  int size = SetLoopDirection(dir);
1091 
1092  for (int id = 0; id < nd/2; id++)
1093  for (int i = 0; i < size; i++)
1094  {
1095  Swap<double>((*this)(id,i), (*this)(nd-1-id,i));
1096  }
1097  kv[dir]->Flip();
1098 }
1099 
1100 void NURBSPatch::SwapDirections(int dir1, int dir2)
1101 {
1102  if (abs(dir1-dir2) == 2)
1103  {
1104  mfem_error("NURBSPatch::SwapDirections :"
1105  " directions 0 and 2 are not supported!");
1106  }
1107 
1108  Array<const KnotVector *> nkv(kv);
1109 
1110  Swap<const KnotVector *>(nkv[dir1], nkv[dir2]);
1111  NURBSPatch *newpatch = new NURBSPatch(nkv, Dim);
1112 
1113  int size = SetLoopDirection(dir1);
1114  newpatch->SetLoopDirection(dir2);
1115 
1116  for (int id = 0; id < nd; id++)
1117  for (int i = 0; i < size; i++)
1118  {
1119  (*newpatch)(id,i) = (*this)(id,i);
1120  }
1121 
1122  swap(newpatch);
1123 }
1124 
1125 void NURBSPatch::Get3DRotationMatrix(double n[], double angle, double r,
1126  DenseMatrix &T)
1127 {
1128  double c, s, c1;
1129  double l2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
1130  double l = sqrt(l2);
1131 
1132  if (fabs(angle) == M_PI_2)
1133  {
1134  s = r*copysign(1., angle);
1135  c = 0.;
1136  c1 = -1.;
1137  }
1138  else if (fabs(angle) == M_PI)
1139  {
1140  s = 0.;
1141  c = -r;
1142  c1 = c - 1.;
1143  }
1144  else
1145  {
1146  s = r*sin(angle);
1147  c = r*cos(angle);
1148  c1 = c - 1.;
1149  }
1150 
1151  T.SetSize(3);
1152 
1153  T(0,0) = (n[0]*n[0] + (n[1]*n[1] + n[2]*n[2])*c)/l2;
1154  T(0,1) = -(n[0]*n[1]*c1)/l2 - (n[2]*s)/l;
1155  T(0,2) = -(n[0]*n[2]*c1)/l2 + (n[1]*s)/l;
1156  T(1,0) = -(n[0]*n[1]*c1)/l2 + (n[2]*s)/l;
1157  T(1,1) = (n[1]*n[1] + (n[0]*n[0] + n[2]*n[2])*c)/l2;
1158  T(1,2) = -(n[1]*n[2]*c1)/l2 - (n[0]*s)/l;
1159  T(2,0) = -(n[0]*n[2]*c1)/l2 - (n[1]*s)/l;
1160  T(2,1) = -(n[1]*n[2]*c1)/l2 + (n[0]*s)/l;
1161  T(2,2) = (n[2]*n[2] + (n[0]*n[0] + n[1]*n[1])*c)/l2;
1162 }
1163 
1164 void NURBSPatch::Rotate3D(double n[], double angle)
1165 {
1166  if (Dim != 4)
1167  {
1168  mfem_error("NURBSPatch::Rotate3D : not a NURBSPatch in 3D!");
1169  }
1170 
1171  DenseMatrix T(3);
1172  Vector x(3), y(NULL, 3);
1173 
1174  Get3DRotationMatrix(n, angle, 1., T);
1175 
1176  int size = 1;
1177  for (int i = 0; i < kv.Size(); i++)
1178  {
1179  size *= kv[i]->GetNCP();
1180  }
1181 
1182  for (int i = 0; i < size; i++)
1183  {
1184  y.SetData(data + i*Dim);
1185  x = y;
1186  T.Mult(x, y);
1187  }
1188 }
1189 
1190 int NURBSPatch::MakeUniformDegree(int degree)
1191 {
1192  int maxd = degree;
1193 
1194  if (maxd == -1)
1195  {
1196  for (int dir = 0; dir < kv.Size(); dir++)
1197  {
1198  maxd = std::max(maxd, kv[dir]->GetOrder());
1199  }
1200  }
1201 
1202  for (int dir = 0; dir < kv.Size(); dir++)
1203  {
1204  if (maxd > kv[dir]->GetOrder())
1205  {
1206  DegreeElevate(dir, maxd - kv[dir]->GetOrder());
1207  }
1208  }
1209 
1210  return maxd;
1211 }
1212 
1213 NURBSPatch *Interpolate(NURBSPatch &p1, NURBSPatch &p2)
1214 {
1215  if (p1.kv.Size() != p2.kv.Size() || p1.Dim != p2.Dim)
1216  {
1217  mfem_error("Interpolate(NURBSPatch &, NURBSPatch &)");
1218  }
1219 
1220  int size = 1, dim = p1.Dim;
1221  Array<const KnotVector *> kv(p1.kv.Size() + 1);
1222 
1223  for (int i = 0; i < p1.kv.Size(); i++)
1224  {
1225  if (p1.kv[i]->GetOrder() < p2.kv[i]->GetOrder())
1226  {
1227  p1.KnotInsert(i, *p2.kv[i]);
1228  p2.KnotInsert(i, *p1.kv[i]);
1229  }
1230  else
1231  {
1232  p2.KnotInsert(i, *p1.kv[i]);
1233  p1.KnotInsert(i, *p2.kv[i]);
1234  }
1235  kv[i] = p1.kv[i];
1236  size *= kv[i]->GetNCP();
1237  }
1238 
1239  KnotVector &nkv = *(new KnotVector(1, 2));
1240  nkv[0] = nkv[1] = 0.0;
1241  nkv[2] = nkv[3] = 1.0;
1242  nkv.GetElements();
1243  kv.Last() = &nkv;
1244 
1245  NURBSPatch *patch = new NURBSPatch(kv, dim);
1246  delete kv.Last();
1247 
1248  for (int i = 0; i < size; i++)
1249  {
1250  for (int d = 0; d < dim; d++)
1251  {
1252  patch->data[i*dim+d] = p1.data[i*dim+d];
1253  patch->data[(i+size)*dim+d] = p2.data[i*dim+d];
1254  }
1255  }
1256 
1257  return patch;
1258 }
1259 
1260 NURBSPatch *Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
1261 {
1262  if (patch.Dim != 4)
1263  {
1264  mfem_error("Revolve3D(NURBSPatch &, double [], double)");
1265  }
1266 
1267  int size = 1, ns;
1268  Array<const KnotVector *> nkv(patch.kv.Size() + 1);
1269 
1270  for (int i = 0; i < patch.kv.Size(); i++)
1271  {
1272  nkv[i] = patch.kv[i];
1273  size *= nkv[i]->GetNCP();
1274  }
1275  ns = 2*times + 1;
1276  KnotVector &lkv = *(new KnotVector(2, ns));
1277  nkv.Last() = &lkv;
1278  lkv[0] = lkv[1] = lkv[2] = 0.0;
1279  for (int i = 1; i < times; i++)
1280  {
1281  lkv[2*i+1] = lkv[2*i+2] = i;
1282  }
1283  lkv[ns] = lkv[ns+1] = lkv[ns+2] = times;
1284  lkv.GetElements();
1285  NURBSPatch *newpatch = new NURBSPatch(nkv, 4);
1286  delete nkv.Last();
1287 
1288  DenseMatrix T(3), T2(3);
1289  Vector u(NULL, 3), v(NULL, 3);
1290 
1291  NURBSPatch::Get3DRotationMatrix(n, ang, 1., T);
1292  double c = cos(ang/2);
1293  NURBSPatch::Get3DRotationMatrix(n, ang/2, 1./c, T2);
1294  T2 *= c;
1295 
1296  double *op = patch.data, *np;
1297  for (int i = 0; i < size; i++)
1298  {
1299  np = newpatch->data + 4*i;
1300  for (int j = 0; j < 4; j++)
1301  {
1302  np[j] = op[j];
1303  }
1304  for (int j = 0; j < times; j++)
1305  {
1306  u.SetData(np);
1307  v.SetData(np += 4*size);
1308  T2.Mult(u, v);
1309  v[3] = c*u[3];
1310  v.SetData(np += 4*size);
1311  T.Mult(u, v);
1312  v[3] = u[3];
1313  }
1314  op += 4;
1315  }
1316 
1317  return newpatch;
1318 }
1319 
1320 
1321 NURBSExtension::NURBSExtension(const NURBSExtension &orig)
1322  : mOrder(orig.mOrder), mOrders(orig.mOrders),
1323  NumOfKnotVectors(orig.NumOfKnotVectors),
1324  NumOfVertices(orig.NumOfVertices),
1325  NumOfElements(orig.NumOfElements),
1326  NumOfBdrElements(orig.NumOfBdrElements),
1327  NumOfDofs(orig.NumOfDofs),
1328  NumOfActiveVertices(orig.NumOfActiveVertices),
1329  NumOfActiveElems(orig.NumOfActiveElems),
1330  NumOfActiveBdrElems(orig.NumOfActiveBdrElems),
1331  NumOfActiveDofs(orig.NumOfActiveDofs),
1332  activeVert(orig.activeVert),
1333  activeElem(orig.activeElem),
1334  activeBdrElem(orig.activeBdrElem),
1335  activeDof(orig.activeDof),
1336  patchTopo(new Mesh(*orig.patchTopo)),
1337  own_topo(true),
1338  edge_to_knot(orig.edge_to_knot),
1339  knotVectors(orig.knotVectors.Size()), // knotVectors are copied in the body
1340  weights(orig.weights),
1341  d_to_d(orig.d_to_d),
1342  master(orig.master),
1343  slave(orig.slave),
1344  v_meshOffsets(orig.v_meshOffsets),
1345  e_meshOffsets(orig.e_meshOffsets),
1346  f_meshOffsets(orig.f_meshOffsets),
1347  p_meshOffsets(orig.p_meshOffsets),
1348  v_spaceOffsets(orig.v_spaceOffsets),
1349  e_spaceOffsets(orig.e_spaceOffsets),
1350  f_spaceOffsets(orig.f_spaceOffsets),
1351  p_spaceOffsets(orig.p_spaceOffsets),
1352  el_dof(orig.el_dof ? new Table(*orig.el_dof) : NULL),
1353  bel_dof(orig.bel_dof ? new Table(*orig.bel_dof) : NULL),
1354  el_to_patch(orig.el_to_patch),
1355  bel_to_patch(orig.bel_to_patch),
1356  el_to_IJK(orig.el_to_IJK),
1357  bel_to_IJK(orig.bel_to_IJK),
1358  patches(orig.patches.Size()) // patches are copied in the body
1359 {
1360  // Copy the knot vectors:
1361  for (int i = 0; i < knotVectors.Size(); i++)
1362  {
1363  knotVectors[i] = new KnotVector(*orig.knotVectors[i]);
1364  }
1365 
1366  // Copy the patches:
1367  for (int p = 0; p < patches.Size(); p++)
1368  {
1369  patches[p] = new NURBSPatch(*orig.patches[p]);
1370  }
1371 }
1372 
1373 NURBSExtension::NURBSExtension(std::istream &input)
1374 {
1375  // Read topology
1376  patchTopo = new Mesh;
1377  patchTopo->LoadPatchTopo(input, edge_to_knot);
1378  own_topo = 1;
1379 
1380  CheckPatches();
1381  // CheckBdrPatches();
1382 
1383  skip_comment_lines(input, '#');
1384 
1385  // Read knotvectors or patches
1386  string ident;
1387  input >> ws >> ident; // 'knotvectors' or 'patches'
1388  if (ident == "knotvectors")
1389  {
1390  input >> NumOfKnotVectors;
1391  knotVectors.SetSize(NumOfKnotVectors);
1392  for (int i = 0; i < NumOfKnotVectors; i++)
1393  {
1394  knotVectors[i] = new KnotVector(input);
1395  }
1396  }
1397  else if (ident == "patches")
1398  {
1399  patches.SetSize(GetNP());
1400  for (int p = 0; p < patches.Size(); p++)
1401  {
1402  skip_comment_lines(input, '#');
1403  patches[p] = new NURBSPatch(input);
1404  }
1405 
1406  NumOfKnotVectors = 0;
1407  for (int i = 0; i < patchTopo->GetNEdges(); i++)
1408  if (NumOfKnotVectors < KnotInd(i))
1409  {
1410  NumOfKnotVectors = KnotInd(i);
1411  }
1412  NumOfKnotVectors++;
1413  knotVectors.SetSize(NumOfKnotVectors);
1414  knotVectors = NULL;
1415 
1416  Array<int> edges, oedge;
1417  for (int p = 0; p < patches.Size(); p++)
1418  {
1419  patchTopo->GetElementEdges(p, edges, oedge);
1420  if (Dimension() == 2)
1421  {
1422  if (knotVectors[KnotInd(edges[0])] == NULL)
1423  {
1424  knotVectors[KnotInd(edges[0])] =
1425  new KnotVector(*patches[p]->GetKV(0));
1426  }
1427  if (knotVectors[KnotInd(edges[1])] == NULL)
1428  {
1429  knotVectors[KnotInd(edges[1])] =
1430  new KnotVector(*patches[p]->GetKV(1));
1431  }
1432  }
1433  else
1434  {
1435  if (knotVectors[KnotInd(edges[0])] == NULL)
1436  {
1437  knotVectors[KnotInd(edges[0])] =
1438  new KnotVector(*patches[p]->GetKV(0));
1439  }
1440  if (knotVectors[KnotInd(edges[3])] == NULL)
1441  {
1442  knotVectors[KnotInd(edges[3])] =
1443  new KnotVector(*patches[p]->GetKV(1));
1444  }
1445  if (knotVectors[KnotInd(edges[8])] == NULL)
1446  {
1447  knotVectors[KnotInd(edges[8])] =
1448  new KnotVector(*patches[p]->GetKV(2));
1449  }
1450  }
1451  }
1452  }
1453  else
1454  {
1455  MFEM_ABORT("invalid section: " << ident);
1456  }
1457 
1458  SetOrdersFromKnotVectors();
1459 
1460  GenerateOffsets();
1461  CountElements();
1462  CountBdrElements();
1463  // NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs
1464 
1465  skip_comment_lines(input, '#');
1466 
1467  // Check for a list of mesh elements
1468  if (patches.Size() == 0)
1469  {
1470  input >> ws >> ident;
1471  }
1472  if (patches.Size() == 0 && ident == "mesh_elements")
1473  {
1474  input >> NumOfActiveElems;
1475  activeElem.SetSize(GetGNE());
1476  activeElem = false;
1477  int glob_elem;
1478  for (int i = 0; i < NumOfActiveElems; i++)
1479  {
1480  input >> glob_elem;
1481  activeElem[glob_elem] = true;
1482  }
1483 
1484  skip_comment_lines(input, '#');
1485  input >> ws >> ident;
1486  }
1487  else
1488  {
1489  NumOfActiveElems = NumOfElements;
1490  activeElem.SetSize(NumOfElements);
1491  activeElem = true;
1492  }
1493 
1494  GenerateActiveVertices();
1495  InitDofMap();
1496  GenerateElementDofTable();
1497  GenerateActiveBdrElems();
1498  GenerateBdrElementDofTable();
1499 
1500  // periodic
1501  if (ident == "periodic")
1502  {
1503  master.Load(input);
1504  slave.Load(input);
1505 
1506  skip_comment_lines(input, '#');
1507  input >> ws >> ident;
1508  }
1509 
1510  if (patches.Size() == 0)
1511  {
1512  // weights
1513  if (ident == "weights")
1514  {
1515  weights.Load(input, GetNDof());
1516  }
1517  else // e.g. ident = "unitweights" or "autoweights"
1518  {
1519  weights.SetSize(GetNDof());
1520  weights = 1.0;
1521  }
1522  }
1523 
1524  // periodic
1525  ConnectBoundaries();
1526 }
1527 
1528 NURBSExtension::NURBSExtension(NURBSExtension *parent, int newOrder)
1529 {
1530  patchTopo = parent->patchTopo;
1531  own_topo = 0;
1532 
1533  parent->edge_to_knot.Copy(edge_to_knot);
1534 
1535  NumOfKnotVectors = parent->GetNKV();
1536  knotVectors.SetSize(NumOfKnotVectors);
1537  const Array<int> &pOrders = parent->GetOrders();
1538  for (int i = 0; i < NumOfKnotVectors; i++)
1539  {
1540  if (newOrder > pOrders[i])
1541  {
1542  knotVectors[i] =
1543  parent->GetKnotVector(i)->DegreeElevate(newOrder - pOrders[i]);
1544  }
1545  else
1546  {
1547  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1548  }
1549  }
1550 
1551  // copy some data from parent
1552  NumOfElements = parent->NumOfElements;
1553  NumOfBdrElements = parent->NumOfBdrElements;
1554 
1555  SetOrdersFromKnotVectors();
1556 
1557  GenerateOffsets(); // dof offsets will be different from parent
1558 
1559  NumOfActiveVertices = parent->NumOfActiveVertices;
1560  NumOfActiveElems = parent->NumOfActiveElems;
1561  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1562  parent->activeVert.Copy(activeVert);
1563  InitDofMap();
1564  parent->activeElem.Copy(activeElem);
1565  parent->activeBdrElem.Copy(activeBdrElem);
1566 
1567  GenerateElementDofTable();
1568  GenerateBdrElementDofTable();
1569 
1570  weights.SetSize(GetNDof());
1571  weights = 1.0;
1572 
1573  // periodic
1574  parent->master.Copy(master);
1575  parent->slave.Copy(slave);
1576  ConnectBoundaries();
1577 }
1578 
1579 NURBSExtension::NURBSExtension(NURBSExtension *parent,
1580  const Array<int> &newOrders)
1581 {
1582  newOrders.Copy(mOrders);
1583  SetOrderFromOrders();
1584 
1585  patchTopo = parent->patchTopo;
1586  own_topo = 0;
1587 
1588  parent->edge_to_knot.Copy(edge_to_knot);
1589 
1590  NumOfKnotVectors = parent->GetNKV();
1591  MFEM_VERIFY(mOrders.Size() == NumOfKnotVectors, "invalid newOrders array");
1592  knotVectors.SetSize(NumOfKnotVectors);
1593  const Array<int> &pOrders = parent->GetOrders();
1594 
1595  for (int i = 0; i < NumOfKnotVectors; i++)
1596  {
1597  if (mOrders[i] > pOrders[i])
1598  {
1599  knotVectors[i] =
1600  parent->GetKnotVector(i)->DegreeElevate(mOrders[i] - pOrders[i]);
1601  }
1602  else
1603  {
1604  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1605  }
1606  }
1607 
1608  // copy some data from parent
1609  NumOfElements = parent->NumOfElements;
1610  NumOfBdrElements = parent->NumOfBdrElements;
1611 
1612  GenerateOffsets(); // dof offsets will be different from parent
1613 
1614  NumOfActiveVertices = parent->NumOfActiveVertices;
1615  NumOfActiveElems = parent->NumOfActiveElems;
1616  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1617  parent->activeVert.Copy(activeVert);
1618  InitDofMap();
1619  parent->activeElem.Copy(activeElem);
1620  parent->activeBdrElem.Copy(activeBdrElem);
1621 
1622  GenerateElementDofTable();
1623  GenerateBdrElementDofTable();
1624 
1625  weights.SetSize(GetNDof());
1626  weights = 1.0;
1627 
1628  parent->master.Copy(master);
1629  parent->slave.Copy(slave);
1630  ConnectBoundaries();
1631 }
1632 
1633 NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
1634 {
1635  NURBSExtension *parent = mesh_array[0]->NURBSext;
1636 
1637  if (!parent->own_topo)
1638  {
1639  mfem_error("NURBSExtension::NURBSExtension :\n"
1640  " parent does not own the patch topology!");
1641  }
1642  patchTopo = parent->patchTopo;
1643  own_topo = 1;
1644  parent->own_topo = 0;
1645 
1646  parent->edge_to_knot.Copy(edge_to_knot);
1647 
1648  parent->GetOrders().Copy(mOrders);
1649  mOrder = parent->GetOrder();
1650 
1651  NumOfKnotVectors = parent->GetNKV();
1652  knotVectors.SetSize(NumOfKnotVectors);
1653  for (int i = 0; i < NumOfKnotVectors; i++)
1654  {
1655  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1656  }
1657 
1658  GenerateOffsets();
1659  CountElements();
1660  CountBdrElements();
1661 
1662  // assuming the meshes define a partitioning of all the elements
1663  NumOfActiveElems = NumOfElements;
1664  activeElem.SetSize(NumOfElements);
1665  activeElem = true;
1666 
1667  GenerateActiveVertices();
1668  InitDofMap();
1669  GenerateElementDofTable();
1670  GenerateActiveBdrElems();
1671  GenerateBdrElementDofTable();
1672 
1673  weights.SetSize(GetNDof());
1674  MergeWeights(mesh_array, num_pieces);
1675 }
1676 
1677 NURBSExtension::~NURBSExtension()
1678 {
1679  if (patches.Size() == 0)
1680  {
1681  delete bel_dof;
1682  delete el_dof;
1683  }
1684 
1685  for (int i = 0; i < knotVectors.Size(); i++)
1686  {
1687  delete knotVectors[i];
1688  }
1689 
1690  for (int i = 0; i < patches.Size(); i++)
1691  {
1692  delete patches[i];
1693  }
1694 
1695  if (own_topo)
1696  {
1697  delete patchTopo;
1698  }
1699 }
1700 
1701 void NURBSExtension::Print(std::ostream &os) const
1702 {
1703  patchTopo->PrintTopo(os, edge_to_knot);
1704  if (patches.Size() == 0)
1705  {
1706  os << "\nknotvectors\n" << NumOfKnotVectors << '\n';
1707  for (int i = 0; i < NumOfKnotVectors; i++)
1708  {
1709  knotVectors[i]->Print(os);
1710  }
1711 
1712  if (NumOfActiveElems < NumOfElements)
1713  {
1714  os << "\nmesh_elements\n" << NumOfActiveElems << '\n';
1715  for (int i = 0; i < NumOfElements; i++)
1716  if (activeElem[i])
1717  {
1718  os << i << '\n';
1719  }
1720  }
1721 
1722  os << "\nweights\n";
1723  weights.Print(os, 1);
1724  }
1725  else
1726  {
1727  os << "\npatches\n";
1728  for (int p = 0; p < patches.Size(); p++)
1729  {
1730  os << "\n# patch " << p << "\n\n";
1731  patches[p]->Print(os);
1732  }
1733  }
1734 }
1735 
1736 void NURBSExtension::PrintCharacteristics(std::ostream &os) const
1737 {
1738  os <<
1739  "NURBS Mesh entity sizes:\n"
1740  "Dimension = " << Dimension() << "\n"
1741  "Unique Orders = ";
1742  Array<int> unique_orders(mOrders);
1743  unique_orders.Sort();
1744  unique_orders.Unique();
1745  unique_orders.Print(os, unique_orders.Size());
1746  os <<
1747  "NumOfKnotVectors = " << GetNKV() << "\n"
1748  "NumOfPatches = " << GetNP() << "\n"
1749  "NumOfBdrPatches = " << GetNBP() << "\n"
1750  "NumOfVertices = " << GetGNV() << "\n"
1751  "NumOfElements = " << GetGNE() << "\n"
1752  "NumOfBdrElements = " << GetGNBE() << "\n"
1753  "NumOfDofs = " << GetNTotalDof() << "\n"
1754  "NumOfActiveVertices = " << GetNV() << "\n"
1755  "NumOfActiveElems = " << GetNE() << "\n"
1756  "NumOfActiveBdrElems = " << GetNBE() << "\n"
1757  "NumOfActiveDofs = " << GetNDof() << '\n';
1758  for (int i = 0; i < NumOfKnotVectors; i++)
1759  {
1760  os << ' ' << i + 1 << ") ";
1761  knotVectors[i]->Print(os);
1762  }
1763  os << endl;
1764 }
1765 
1766 void NURBSExtension::PrintFunctions(const char *basename, int samples) const
1767 {
1768  std::ofstream os;
1769  for (int i = 0; i < NumOfKnotVectors; i++)
1770  {
1771  std::ostringstream filename;
1772  filename << basename<<"_"<<i<<".dat";
1773  os.open(filename.str().c_str());
1774  knotVectors[i]->PrintFunctions(os,samples);
1775  os.close();
1776  }
1777 }
1778 
1779 void NURBSExtension::InitDofMap()
1780 {
1781  master.SetSize(0);
1782  slave.SetSize(0);
1783  d_to_d.SetSize(0);
1784 }
1785 
1786 void NURBSExtension::ConnectBoundaries(Array<int> &bnds0, Array<int> &bnds1)
1787 {
1788  bnds0.Copy(master);
1789  bnds1.Copy(slave);
1790  ConnectBoundaries();
1791 }
1792 
1793 void NURBSExtension::ConnectBoundaries()
1794 {
1795  if (master.Size() != slave.Size())
1796  {
1797  mfem_error("NURBSExtension::ConnectBoundaries() boundary lists not of equal size");
1798  }
1799  if (master.Size() == 0 ) { return; }
1800 
1801  // Connect
1802  d_to_d.SetSize(NumOfDofs);
1803  for (int i = 0; i < NumOfDofs; i++) { d_to_d[i] = i; }
1804 
1805  for (int i = 0; i < master.Size(); i++)
1806  {
1807  if (Dimension() == 2)
1808  {
1809  ConnectBoundaries2D(master[i], slave[i]);
1810  }
1811  else
1812  {
1813  ConnectBoundaries3D(master[i], slave[i]);
1814  }
1815  }
1816 
1817  // Finalize
1818  if (el_dof) { delete el_dof; }
1819  if (bel_dof) { delete bel_dof; }
1820  GenerateElementDofTable();
1821  GenerateBdrElementDofTable();
1822 }
1823 
1824 void NURBSExtension::ConnectBoundaries2D(int bnd0, int bnd1)
1825 {
1826  int idx0 = -1, idx1 = -1;
1827  for (int b = 0; b < GetNBP(); b++)
1828  {
1829  if (bnd0 == patchTopo->GetBdrAttribute(b)) { idx0 = b; }
1830  if (bnd1 == patchTopo->GetBdrAttribute(b)) { idx1 = b; }
1831  }
1832  MFEM_VERIFY(idx0 != -1,"Bdr 0 not found");
1833  MFEM_VERIFY(idx1 != -1,"Bdr 1 not found");
1834 
1835  NURBSPatchMap p2g0(this);
1836  NURBSPatchMap p2g1(this);
1837 
1838  int okv0[1],okv1[1];
1839  const KnotVector *kv0[1],*kv1[1];
1840 
1841  p2g0.SetBdrPatchDofMap(idx0, kv0, okv0);
1842  p2g1.SetBdrPatchDofMap(idx1, kv1, okv1);
1843 
1844  int nx = p2g0.nx();
1845  int nks0 = kv0[0]->GetNKS();
1846 
1847 #ifdef MFEM_DEBUG
1848  bool compatible = true;
1849  if (p2g0.nx() != p2g1.nx()) { compatible = false; }
1850  if (kv0[0]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
1851  if (kv0[0]->GetOrder() != kv1[0]->GetOrder()) { compatible = false; }
1852 
1853  if (!compatible)
1854  {
1855  mfem::out<<p2g0.nx()<<" "<<p2g1.nx()<<endl;
1856  mfem::out<<kv0[0]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
1857  mfem::out<<kv0[0]->GetOrder()<<" "<<kv1[0]->GetOrder()<<endl;
1858  mfem_error("NURBS boundaries not compatible");
1859  }
1860 #endif
1861 
1862  for (int i = 0; i < nks0; i++)
1863  {
1864  if (kv0[0]->isElement(i))
1865  {
1866  if (!kv1[0]->isElement(i)) { mfem_error("isElement does not match"); }
1867  for (int ii = 0; ii <= kv0[0]->GetOrder(); ii++)
1868  {
1869  int ii0 = (okv0[0] >= 0) ? (i+ii) : (nx-i-ii);
1870  int ii1 = (okv1[0] >= 0) ? (i+ii) : (nx-i-ii);
1871 
1872  d_to_d[p2g0(ii0)] = d_to_d[p2g1(ii1)];
1873  }
1874 
1875  }
1876  }
1877 
1878  // Clean d_to_d
1879  Array<int> tmp(d_to_d.Size()+1);
1880  tmp = 0;
1881 
1882  for (int i = 0; i < d_to_d.Size(); i++)
1883  {
1884  tmp[d_to_d[i]] = 1;
1885  }
1886 
1887  int cnt = 0;
1888  for (int i = 0; i < tmp.Size(); i++)
1889  {
1890  if (tmp[i] == 1) { tmp[i] = cnt++; }
1891  }
1892  NumOfDofs = cnt;
1893 
1894  for (int i = 0; i < d_to_d.Size(); i++)
1895  {
1896  d_to_d[i] = tmp[d_to_d[i]];
1897  }
1898 }
1899 
1900 void NURBSExtension::ConnectBoundaries3D(int bnd0, int bnd1)
1901 {
1902  int idx0 = -1, idx1 = -1;
1903  for (int b = 0; b < GetNBP(); b++)
1904  {
1905  if (bnd0 == patchTopo->GetBdrAttribute(b)) { idx0 = b; }
1906  if (bnd1 == patchTopo->GetBdrAttribute(b)) { idx1 = b; }
1907  }
1908  MFEM_VERIFY(idx0 != -1,"Bdr 0 not found");
1909  MFEM_VERIFY(idx1 != -1,"Bdr 1 not found");
1910 
1911  NURBSPatchMap p2g0(this);
1912  NURBSPatchMap p2g1(this);
1913 
1914  int okv0[2],okv1[2];
1915  const KnotVector *kv0[2],*kv1[2];
1916 
1917  p2g0.SetBdrPatchDofMap(idx0, kv0, okv0);
1918  p2g1.SetBdrPatchDofMap(idx1, kv1, okv1);
1919 
1920  int nx = p2g0.nx();
1921  int ny = p2g0.ny();
1922 
1923  int nks0 = kv0[0]->GetNKS();
1924  int nks1 = kv0[1]->GetNKS();
1925 
1926 #ifdef MFEM_DEBUG
1927  bool compatible = true;
1928  if (p2g0.nx() != p2g1.nx()) { compatible = false; }
1929  if (p2g0.ny() != p2g1.ny()) { compatible = false; }
1930 
1931  if (kv0[0]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
1932  if (kv0[1]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
1933 
1934  if (kv0[0]->GetOrder() != kv1[0]->GetOrder()) { compatible = false; }
1935  if (kv0[1]->GetOrder() != kv1[1]->GetOrder()) { compatible = false; }
1936 
1937  if (!compatible)
1938  {
1939  mfem::out<<p2g0.nx()<<" "<<p2g1.nx()<<endl;
1940  mfem::out<<p2g0.ny()<<" "<<p2g1.ny()<<endl;
1941 
1942  mfem::out<<kv0[0]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
1943  mfem::out<<kv0[1]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
1944 
1945  mfem::out<<kv0[0]->GetOrder()<<" "<<kv1[0]->GetOrder()<<endl;
1946  mfem::out<<kv0[1]->GetOrder()<<" "<<kv1[1]->GetOrder()<<endl;
1947  mfem_error("NURBS boundaries not compatible");
1948  }
1949 #endif
1950 
1951  for (int j = 0; j < nks1; j++)
1952  {
1953  if (kv0[1]->isElement(j))
1954  {
1955  if (!kv1[1]->isElement(j)) { mfem_error("isElement does not match #1"); }
1956  for (int i = 0; i < nks0; i++)
1957  {
1958  if (kv0[0]->isElement(i))
1959  {
1960  if (!kv1[0]->isElement(i)) { mfem_error("isElement does not match #0"); }
1961  for (int jj = 0; jj <= kv0[1]->GetOrder(); jj++)
1962  {
1963  int jj0 = (okv0[1] >= 0) ? (j+jj) : (ny-j-jj);
1964  int jj1 = (okv1[1] >= 0) ? (j+jj) : (ny-j-jj);
1965 
1966  for (int ii = 0; ii <= kv0[0]->GetOrder(); ii++)
1967  {
1968  int ii0 = (okv0[0] >= 0) ? (i+ii) : (nx-i-ii);
1969  int ii1 = (okv1[0] >= 0) ? (i+ii) : (nx-i-ii);
1970 
1971  d_to_d[p2g0(ii0,jj0)] = d_to_d[p2g1(ii1,jj1)];
1972  }
1973  }
1974  }
1975  }
1976  }
1977  }
1978 
1979  // Clean d_to_d
1980  Array<int> tmp(d_to_d.Size()+1);
1981  tmp = 0;
1982 
1983  for (int i = 0; i < d_to_d.Size(); i++)
1984  {
1985  tmp[d_to_d[i]] = 1;
1986  }
1987 
1988  int cnt = 0;
1989  for (int i = 0; i < tmp.Size(); i++)
1990  {
1991  if (tmp[i] == 1) { tmp[i] = cnt++; }
1992  }
1993  NumOfDofs = cnt;
1994 
1995  for (int i = 0; i < d_to_d.Size(); i++)
1996  {
1997  d_to_d[i] = tmp[d_to_d[i]];
1998  }
1999 }
2000 
2002 {
2003  int vert[8], nv, g_el, nx, ny, nz, dim = Dimension();
2004 
2005  NURBSPatchMap p2g(this);
2006  const KnotVector *kv[3];
2007 
2008  g_el = 0;
2009  activeVert.SetSize(GetGNV());
2010  activeVert = -1;
2011  for (int p = 0; p < GetNP(); p++)
2012  {
2013  p2g.SetPatchVertexMap(p, kv);
2014 
2015  nx = p2g.nx();
2016  ny = p2g.ny();
2017  nz = (dim == 3) ? p2g.nz() : 1;
2018 
2019  for (int k = 0; k < nz; k++)
2020  {
2021  for (int j = 0; j < ny; j++)
2022  {
2023  for (int i = 0; i < nx; i++)
2024  {
2025  if (activeElem[g_el])
2026  {
2027  if (dim == 2)
2028  {
2029  vert[0] = p2g(i, j );
2030  vert[1] = p2g(i+1,j );
2031  vert[2] = p2g(i+1,j+1);
2032  vert[3] = p2g(i, j+1);
2033  nv = 4;
2034  }
2035  else
2036  {
2037  vert[0] = p2g(i, j, k);
2038  vert[1] = p2g(i+1,j, k);
2039  vert[2] = p2g(i+1,j+1,k);
2040  vert[3] = p2g(i, j+1,k);
2041 
2042  vert[4] = p2g(i, j, k+1);
2043  vert[5] = p2g(i+1,j, k+1);
2044  vert[6] = p2g(i+1,j+1,k+1);
2045  vert[7] = p2g(i, j+1,k+1);
2046  nv = 8;
2047  }
2048 
2049  for (int v = 0; v < nv; v++)
2050  {
2051  activeVert[vert[v]] = 1;
2052  }
2053  }
2054  g_el++;
2055  }
2056  }
2057  }
2058  }
2059 
2060  NumOfActiveVertices = 0;
2061  for (int i = 0; i < GetGNV(); i++)
2062  if (activeVert[i] == 1)
2063  {
2064  activeVert[i] = NumOfActiveVertices++;
2065  }
2066 }
2067 
2069 {
2070  int dim = Dimension();
2071  Array<KnotVector *> kv(dim);
2072 
2073  activeBdrElem.SetSize(GetGNBE());
2074  if (GetGNE() == GetNE())
2075  {
2076  activeBdrElem = true;
2077  NumOfActiveBdrElems = GetGNBE();
2078  return;
2079  }
2080  activeBdrElem = false;
2081  NumOfActiveBdrElems = 0;
2082  // the mesh will generate the actual boundary including boundary
2083  // elements that are not on boundary patches. we use this for
2084  // visualization of processor boundaries
2085 
2086  // TODO: generate actual boundary?
2087 }
2088 
2089 
2090 void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces)
2091 {
2092  Array<int> lelem_elem;
2093 
2094  for (int i = 0; i < num_pieces; i++)
2095  {
2096  NURBSExtension *lext = mesh_array[i]->NURBSext;
2097 
2098  lext->GetElementLocalToGlobal(lelem_elem);
2099 
2100  for (int lel = 0; lel < lext->GetNE(); lel++)
2101  {
2102  int gel = lelem_elem[lel];
2103 
2104  int nd = el_dof->RowSize(gel);
2105  int *gdofs = el_dof->GetRow(gel);
2106  int *ldofs = lext->el_dof->GetRow(lel);
2107  for (int j = 0; j < nd; j++)
2108  {
2109  weights(gdofs[j]) = lext->weights(ldofs[j]);
2110  }
2111  }
2112  }
2113 }
2114 
2116  GridFunction *gf_array[], int num_pieces, GridFunction &merged)
2117 {
2118  FiniteElementSpace *gfes = merged.FESpace();
2119  Array<int> lelem_elem, dofs;
2120  Vector lvec;
2121 
2122  for (int i = 0; i < num_pieces; i++)
2123  {
2124  FiniteElementSpace *lfes = gf_array[i]->FESpace();
2125  NURBSExtension *lext = lfes->GetMesh()->NURBSext;
2126 
2127  lext->GetElementLocalToGlobal(lelem_elem);
2128 
2129  for (int lel = 0; lel < lext->GetNE(); lel++)
2130  {
2131  lfes->GetElementVDofs(lel, dofs);
2132  gf_array[i]->GetSubVector(dofs, lvec);
2133 
2134  gfes->GetElementVDofs(lelem_elem[lel], dofs);
2135  merged.SetSubVector(dofs, lvec);
2136  }
2137  }
2138 }
2139 
2141 {
2142  Array<int> edges;
2143  Array<int> oedge;
2144 
2145  for (int p = 0; p < GetNP(); p++)
2146  {
2147  patchTopo->GetElementEdges(p, edges, oedge);
2148 
2149  for (int i = 0; i < edges.Size(); i++)
2150  {
2151  edges[i] = edge_to_knot[edges[i]];
2152  if (oedge[i] < 0)
2153  {
2154  edges[i] = -1 - edges[i];
2155  }
2156  }
2157 
2158  if ((Dimension() == 2 &&
2159  (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
2160 
2161  (Dimension() == 3 &&
2162  (edges[0] != edges[2] || edges[0] != edges[4] ||
2163  edges[0] != edges[6] || edges[1] != edges[3] ||
2164  edges[1] != edges[5] || edges[1] != edges[7] ||
2165  edges[8] != edges[9] || edges[8] != edges[10] ||
2166  edges[8] != edges[11])))
2167  {
2168  mfem::err << "NURBSExtension::CheckPatch (patch = " << p
2169  << ")\n Inconsistent edge-to-knot mapping!\n";
2170  mfem_error();
2171  }
2172 
2173  if ((Dimension() == 2 &&
2174  (edges[0] < 0 || edges[1] < 0)) ||
2175 
2176  (Dimension() == 3 &&
2177  (edges[0] < 0 || edges[3] < 0 || edges[8] < 0)))
2178  {
2179  mfem::err << "NURBSExtension::CheckPatch (patch = " << p
2180  << ") : Bad orientation!\n";
2181  mfem_error();
2182  }
2183  }
2184 }
2185 
2187 {
2188  Array<int> edges;
2189  Array<int> oedge;
2190 
2191  for (int p = 0; p < GetNBP(); p++)
2192  {
2193  patchTopo->GetBdrElementEdges(p, edges, oedge);
2194 
2195  for (int i = 0; i < edges.Size(); i++)
2196  {
2197  edges[i] = edge_to_knot[edges[i]];
2198  if (oedge[i] < 0)
2199  {
2200  edges[i] = -1 - edges[i];
2201  }
2202  }
2203 
2204  if ((Dimension() == 2 && (edges[0] < 0)) ||
2205  (Dimension() == 3 && (edges[0] < 0 || edges[1] < 0)))
2206  {
2207  mfem::err << "NURBSExtension::CheckBdrPatch (boundary patch = "
2208  << p << ") : Bad orientation!\n";
2209  mfem_error();
2210  }
2211  }
2212 }
2213 
2215 {
2216  Array<int> edges, orient;
2217 
2218  kv.SetSize(Dimension());
2219  patchTopo->GetElementEdges(p, edges, orient);
2220 
2221  if (Dimension() == 2)
2222  {
2223  kv[0] = KnotVec(edges[0]);
2224  kv[1] = KnotVec(edges[1]);
2225  }
2226  else
2227  {
2228  kv[0] = KnotVec(edges[0]);
2229  kv[1] = KnotVec(edges[3]);
2230  kv[2] = KnotVec(edges[8]);
2231  }
2232 }
2233 
2235 const
2236 {
2237  Array<int> edges, orient;
2238 
2239  kv.SetSize(Dimension());
2240  patchTopo->GetElementEdges(p, edges, orient);
2241 
2242  if (Dimension() == 2)
2243  {
2244  kv[0] = KnotVec(edges[0]);
2245  kv[1] = KnotVec(edges[1]);
2246  }
2247  else
2248  {
2249  kv[0] = KnotVec(edges[0]);
2250  kv[1] = KnotVec(edges[3]);
2251  kv[2] = KnotVec(edges[8]);
2252  }
2253 }
2254 
2256 {
2257  Array<int> edges;
2258  Array<int> orient;
2259 
2260  kv.SetSize(Dimension() - 1);
2261  patchTopo->GetBdrElementEdges(p, edges, orient);
2262 
2263  if (Dimension() == 2)
2264  {
2265  kv[0] = KnotVec(edges[0]);
2266  }
2267  else
2268  {
2269  kv[0] = KnotVec(edges[0]);
2270  kv[1] = KnotVec(edges[1]);
2271  }
2272 }
2273 
2275  int p, Array<const KnotVector *> &kv) const
2276 {
2277  Array<int> edges;
2278  Array<int> orient;
2279 
2280  kv.SetSize(Dimension() - 1);
2281  patchTopo->GetBdrElementEdges(p, edges, orient);
2282 
2283  if (Dimension() == 2)
2284  {
2285  kv[0] = KnotVec(edges[0]);
2286  }
2287  else
2288  {
2289  kv[0] = KnotVec(edges[0]);
2290  kv[1] = KnotVec(edges[1]);
2291  }
2292 }
2293 
2295 {
2296  MFEM_VERIFY(mOrders.Size() > 0, "");
2297  mOrder = mOrders[0];
2298  for (int i = 1; i < mOrders.Size(); i++)
2299  {
2300  if (mOrders[i] != mOrder)
2301  {
2303  return;
2304  }
2305  }
2306 }
2307 
2309 {
2310  mOrders.SetSize(NumOfKnotVectors);
2311  for (int i = 0; i < NumOfKnotVectors; i++)
2312  {
2313  mOrders[i] = knotVectors[i]->GetOrder();
2314  }
2315  SetOrderFromOrders();
2316 }
2317 
2319 {
2320  int nv = patchTopo->GetNV();
2321  int ne = patchTopo->GetNEdges();
2322  int nf = patchTopo->GetNFaces();
2323  int np = patchTopo->GetNE();
2324  int meshCounter, spaceCounter, dim = Dimension();
2325 
2326  Array<int> edges;
2327  Array<int> orient;
2328 
2329  v_meshOffsets.SetSize(nv);
2330  e_meshOffsets.SetSize(ne);
2331  f_meshOffsets.SetSize(nf);
2332  p_meshOffsets.SetSize(np);
2333 
2334  v_spaceOffsets.SetSize(nv);
2335  e_spaceOffsets.SetSize(ne);
2336  f_spaceOffsets.SetSize(nf);
2337  p_spaceOffsets.SetSize(np);
2338 
2339  // Get vertex offsets
2340  for (meshCounter = 0; meshCounter < nv; meshCounter++)
2341  {
2342  v_meshOffsets[meshCounter] = meshCounter;
2343  v_spaceOffsets[meshCounter] = meshCounter;
2344  }
2345  spaceCounter = meshCounter;
2346 
2347  // Get edge offsets
2348  for (int e = 0; e < ne; e++)
2349  {
2350  e_meshOffsets[e] = meshCounter;
2351  e_spaceOffsets[e] = spaceCounter;
2352  meshCounter += KnotVec(e)->GetNE() - 1;
2353  spaceCounter += KnotVec(e)->GetNCP() - 2;
2354  }
2355 
2356  // Get face offsets
2357  for (int f = 0; f < nf; f++)
2358  {
2359  f_meshOffsets[f] = meshCounter;
2360  f_spaceOffsets[f] = spaceCounter;
2361 
2362  patchTopo->GetFaceEdges(f, edges, orient);
2363 
2364  meshCounter +=
2365  (KnotVec(edges[0])->GetNE() - 1) *
2366  (KnotVec(edges[1])->GetNE() - 1);
2367  spaceCounter +=
2368  (KnotVec(edges[0])->GetNCP() - 2) *
2369  (KnotVec(edges[1])->GetNCP() - 2);
2370  }
2371 
2372  // Get patch offsets
2373  for (int p = 0; p < np; p++)
2374  {
2375  p_meshOffsets[p] = meshCounter;
2376  p_spaceOffsets[p] = spaceCounter;
2377 
2378  patchTopo->GetElementEdges(p, edges, orient);
2379 
2380  if (dim == 2)
2381  {
2382  meshCounter +=
2383  (KnotVec(edges[0])->GetNE() - 1) *
2384  (KnotVec(edges[1])->GetNE() - 1);
2385  spaceCounter +=
2386  (KnotVec(edges[0])->GetNCP() - 2) *
2387  (KnotVec(edges[1])->GetNCP() - 2);
2388  }
2389  else
2390  {
2391  meshCounter +=
2392  (KnotVec(edges[0])->GetNE() - 1) *
2393  (KnotVec(edges[3])->GetNE() - 1) *
2394  (KnotVec(edges[8])->GetNE() - 1);
2395  spaceCounter +=
2396  (KnotVec(edges[0])->GetNCP() - 2) *
2397  (KnotVec(edges[3])->GetNCP() - 2) *
2398  (KnotVec(edges[8])->GetNCP() - 2);
2399  }
2400  }
2401  NumOfVertices = meshCounter;
2402  NumOfDofs = spaceCounter;
2403 }
2404 
2406 {
2407  int dim = Dimension();
2408  Array<const KnotVector *> kv(dim);
2409 
2410  NumOfElements = 0;
2411  for (int p = 0; p < GetNP(); p++)
2412  {
2413  GetPatchKnotVectors(p, kv);
2414 
2415  int ne = kv[0]->GetNE();
2416  for (int d = 1; d < dim; d++)
2417  {
2418  ne *= kv[d]->GetNE();
2419  }
2420 
2421  NumOfElements += ne;
2422  }
2423 }
2424 
2426 {
2427  int dim = Dimension() - 1;
2428  Array<KnotVector *> kv(dim);
2429 
2430  NumOfBdrElements = 0;
2431  for (int p = 0; p < GetNBP(); p++)
2432  {
2433  GetBdrPatchKnotVectors(p, kv);
2434 
2435  int ne = kv[0]->GetNE();
2436  for (int d = 1; d < dim; d++)
2437  {
2438  ne *= kv[d]->GetNE();
2439  }
2440 
2441  NumOfBdrElements += ne;
2442  }
2443 }
2444 
2446 {
2447  elements.SetSize(GetNE());
2448 
2449  if (Dimension() == 2)
2450  {
2451  Get2DElementTopo(elements);
2452  }
2453  else
2454  {
2455  Get3DElementTopo(elements);
2456  }
2457 }
2458 
2460 {
2461  int el = 0;
2462  int eg = 0;
2463  int ind[4];
2464  NURBSPatchMap p2g(this);
2465  const KnotVector *kv[2];
2466 
2467  for (int p = 0; p < GetNP(); p++)
2468  {
2469  p2g.SetPatchVertexMap(p, kv);
2470  int nx = p2g.nx();
2471  int ny = p2g.ny();
2472 
2473  int patch_attr = patchTopo->GetAttribute(p);
2474 
2475  for (int j = 0; j < ny; j++)
2476  {
2477  for (int i = 0; i < nx; i++)
2478  {
2479  if (activeElem[eg])
2480  {
2481  ind[0] = activeVert[p2g(i, j )];
2482  ind[1] = activeVert[p2g(i+1,j )];
2483  ind[2] = activeVert[p2g(i+1,j+1)];
2484  ind[3] = activeVert[p2g(i, j+1)];
2485 
2486  elements[el] = new Quadrilateral(ind, patch_attr);
2487  el++;
2488  }
2489  eg++;
2490  }
2491  }
2492  }
2493 }
2494 
2496 {
2497  int el = 0;
2498  int eg = 0;
2499  int ind[8];
2500  NURBSPatchMap p2g(this);
2501  const KnotVector *kv[3];
2502 
2503  for (int p = 0; p < GetNP(); p++)
2504  {
2505  p2g.SetPatchVertexMap(p, kv);
2506  int nx = p2g.nx();
2507  int ny = p2g.ny();
2508  int nz = p2g.nz();
2509 
2510  int patch_attr = patchTopo->GetAttribute(p);
2511 
2512  for (int k = 0; k < nz; k++)
2513  {
2514  for (int j = 0; j < ny; j++)
2515  {
2516  for (int i = 0; i < nx; i++)
2517  {
2518  if (activeElem[eg])
2519  {
2520  ind[0] = activeVert[p2g(i, j, k)];
2521  ind[1] = activeVert[p2g(i+1,j, k)];
2522  ind[2] = activeVert[p2g(i+1,j+1,k)];
2523  ind[3] = activeVert[p2g(i, j+1,k)];
2524 
2525  ind[4] = activeVert[p2g(i, j, k+1)];
2526  ind[5] = activeVert[p2g(i+1,j, k+1)];
2527  ind[6] = activeVert[p2g(i+1,j+1,k+1)];
2528  ind[7] = activeVert[p2g(i, j+1,k+1)];
2529 
2530  elements[el] = new Hexahedron(ind, patch_attr);
2531  el++;
2532  }
2533  eg++;
2534  }
2535  }
2536  }
2537  }
2538 }
2539 
2541 {
2542  boundary.SetSize(GetNBE());
2543 
2544  if (Dimension() == 2)
2545  {
2546  Get2DBdrElementTopo(boundary);
2547  }
2548  else
2549  {
2550  Get3DBdrElementTopo(boundary);
2551  }
2552 }
2553 
2555 {
2556  int g_be, l_be;
2557  int ind[2], okv[1];
2558  NURBSPatchMap p2g(this);
2559  const KnotVector *kv[1];
2560 
2561  g_be = l_be = 0;
2562  for (int b = 0; b < GetNBP(); b++)
2563  {
2564  p2g.SetBdrPatchVertexMap(b, kv, okv);
2565  int nx = p2g.nx();
2566 
2567  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
2568 
2569  for (int i = 0; i < nx; i++)
2570  {
2571  if (activeBdrElem[g_be])
2572  {
2573  int i_ = (okv[0] >= 0) ? i : (nx - 1 - i);
2574  ind[0] = activeVert[p2g[i_ ]];
2575  ind[1] = activeVert[p2g[i_+1]];
2576 
2577  boundary[l_be] = new Segment(ind, bdr_patch_attr);
2578  l_be++;
2579  }
2580  g_be++;
2581  }
2582  }
2583 }
2584 
2586 {
2587  int g_be, l_be;
2588  int ind[4], okv[2];
2589  NURBSPatchMap p2g(this);
2590  const KnotVector *kv[2];
2591 
2592  g_be = l_be = 0;
2593  for (int b = 0; b < GetNBP(); b++)
2594  {
2595  p2g.SetBdrPatchVertexMap(b, kv, okv);
2596  int nx = p2g.nx();
2597  int ny = p2g.ny();
2598 
2599  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
2600 
2601  for (int j = 0; j < ny; j++)
2602  {
2603  int j_ = (okv[1] >= 0) ? j : (ny - 1 - j);
2604  for (int i = 0; i < nx; i++)
2605  {
2606  if (activeBdrElem[g_be])
2607  {
2608  int i_ = (okv[0] >= 0) ? i : (nx - 1 - i);
2609  ind[0] = activeVert[p2g(i_, j_ )];
2610  ind[1] = activeVert[p2g(i_+1,j_ )];
2611  ind[2] = activeVert[p2g(i_+1,j_+1)];
2612  ind[3] = activeVert[p2g(i_, j_+1)];
2613 
2614  boundary[l_be] = new Quadrilateral(ind, bdr_patch_attr);
2615  l_be++;
2616  }
2617  g_be++;
2618  }
2619  }
2620  }
2621 }
2622 
2624 {
2625  activeDof.SetSize(GetNTotalDof());
2626  activeDof = 0;
2627 
2628  if (Dimension() == 2)
2629  {
2630  Generate2DElementDofTable();
2631  }
2632  else
2633  {
2634  Generate3DElementDofTable();
2635  }
2636 
2637  NumOfActiveDofs = 0;
2638  for (int d = 0; d < GetNTotalDof(); d++)
2639  if (activeDof[d])
2640  {
2641  NumOfActiveDofs++;
2642  activeDof[d] = NumOfActiveDofs;
2643  }
2644 
2645  int *dof = el_dof->GetJ();
2646  int ndof = el_dof->Size_of_connections();
2647  for (int i = 0; i < ndof; i++)
2648  {
2649  dof[i] = activeDof[dof[i]] - 1;
2650  }
2651 }
2652 
2654 {
2655  int el = 0;
2656  int eg = 0;
2657  const KnotVector *kv[2];
2658  NURBSPatchMap p2g(this);
2659 
2660  Array<Connection> el_dof_list;
2661  el_to_patch.SetSize(NumOfActiveElems);
2662  el_to_IJK.SetSize(NumOfActiveElems, 2);
2663 
2664  for (int p = 0; p < GetNP(); p++)
2665  {
2666  p2g.SetPatchDofMap(p, kv);
2667 
2668  // Load dofs
2669  const int ord0 = kv[0]->GetOrder();
2670  const int ord1 = kv[1]->GetOrder();
2671  for (int j = 0; j < kv[1]->GetNKS(); j++)
2672  {
2673  if (kv[1]->isElement(j))
2674  {
2675  for (int i = 0; i < kv[0]->GetNKS(); i++)
2676  {
2677  if (kv[0]->isElement(i))
2678  {
2679  if (activeElem[eg])
2680  {
2681  Connection conn(el,0);
2682  for (int jj = 0; jj <= ord1; jj++)
2683  {
2684  for (int ii = 0; ii <= ord0; ii++)
2685  {
2686  conn.to = DofMap(p2g(i+ii,j+jj));
2687  activeDof[conn.to] = 1;
2688  el_dof_list.Append(conn);
2689  }
2690  }
2691  el_to_patch[el] = p;
2692  el_to_IJK(el,0) = i;
2693  el_to_IJK(el,1) = j;
2694 
2695  el++;
2696  }
2697  eg++;
2698  }
2699  }
2700  }
2701  }
2702  }
2703  // We must NOT sort el_dof_list in this case.
2704  el_dof = new Table(NumOfActiveElems, el_dof_list);
2705 }
2706 
2708 {
2709  int el = 0;
2710  int eg = 0;
2711  const KnotVector *kv[3];
2712  NURBSPatchMap p2g(this);
2713 
2714  Array<Connection> el_dof_list;
2715  el_to_patch.SetSize(NumOfActiveElems);
2716  el_to_IJK.SetSize(NumOfActiveElems, 3);
2717 
2718  for (int p = 0; p < GetNP(); p++)
2719  {
2720  p2g.SetPatchDofMap(p, kv);
2721 
2722  // Load dofs
2723  const int ord0 = kv[0]->GetOrder();
2724  const int ord1 = kv[1]->GetOrder();
2725  const int ord2 = kv[2]->GetOrder();
2726  for (int k = 0; k < kv[2]->GetNKS(); k++)
2727  {
2728  if (kv[2]->isElement(k))
2729  {
2730  for (int j = 0; j < kv[1]->GetNKS(); j++)
2731  {
2732  if (kv[1]->isElement(j))
2733  {
2734  for (int i = 0; i < kv[0]->GetNKS(); i++)
2735  {
2736  if (kv[0]->isElement(i))
2737  {
2738  if (activeElem[eg])
2739  {
2740  Connection conn(el,0);
2741  for (int kk = 0; kk <= ord2; kk++)
2742  {
2743  for (int jj = 0; jj <= ord1; jj++)
2744  {
2745  for (int ii = 0; ii <= ord0; ii++)
2746  {
2747  conn.to = DofMap(p2g(i+ii, j+jj, k+kk));
2748  activeDof[conn.to] = 1;
2749  el_dof_list.Append(conn);
2750  }
2751  }
2752  }
2753 
2754  el_to_patch[el] = p;
2755  el_to_IJK(el,0) = i;
2756  el_to_IJK(el,1) = j;
2757  el_to_IJK(el,2) = k;
2758 
2759  el++;
2760  }
2761  eg++;
2762  }
2763  }
2764  }
2765  }
2766  }
2767  }
2768  }
2769  // We must NOT sort el_dof_list in this case.
2770  el_dof = new Table(NumOfActiveElems, el_dof_list);
2771 }
2772 
2774 {
2775  if (Dimension() == 2)
2776  {
2777  Generate2DBdrElementDofTable();
2778  }
2779  else
2780  {
2781  Generate3DBdrElementDofTable();
2782  }
2783 
2784  int *dof = bel_dof->GetJ();
2785  int ndof = bel_dof->Size_of_connections();
2786  for (int i = 0; i < ndof; i++)
2787  {
2788  dof[i] = activeDof[dof[i]] - 1;
2789  }
2790 }
2791 
2793 {
2794  int gbe = 0;
2795  int lbe = 0, okv[1];
2796  const KnotVector *kv[1];
2797  NURBSPatchMap p2g(this);
2798 
2799  Array<Connection> bel_dof_list;
2800  bel_to_patch.SetSize(NumOfActiveBdrElems);
2801  bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
2802 
2803  for (int b = 0; b < GetNBP(); b++)
2804  {
2805  p2g.SetBdrPatchDofMap(b, kv, okv);
2806  const int nx = p2g.nx(); // NCP-1
2807 
2808  // Load dofs
2809  const int nks0 = kv[0]->GetNKS();
2810  const int ord0 = kv[0]->GetOrder();
2811  for (int i = 0; i < nks0; i++)
2812  {
2813  if (kv[0]->isElement(i))
2814  {
2815  if (activeBdrElem[gbe])
2816  {
2817  Connection conn(lbe,0);
2818  for (int ii = 0; ii <= ord0; ii++)
2819  {
2820  conn.to = DofMap(p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)]);
2821  bel_dof_list.Append(conn);
2822  }
2823  bel_to_patch[lbe] = b;
2824  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2825  lbe++;
2826  }
2827  gbe++;
2828  }
2829  }
2830  }
2831  // We must NOT sort bel_dof_list in this case.
2832  bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
2833 }
2834 
2836 {
2837  int gbe = 0;
2838  int lbe = 0, okv[2];
2839  const KnotVector *kv[2];
2840  NURBSPatchMap p2g(this);
2841 
2842  Array<Connection> bel_dof_list;
2843  bel_to_patch.SetSize(NumOfActiveBdrElems);
2844  bel_to_IJK.SetSize(NumOfActiveBdrElems, 2);
2845 
2846  for (int b = 0; b < GetNBP(); b++)
2847  {
2848  p2g.SetBdrPatchDofMap(b, kv, okv);
2849  const int nx = p2g.nx(); // NCP0-1
2850  const int ny = p2g.ny(); // NCP1-1
2851 
2852  // Load dofs
2853  const int nks0 = kv[0]->GetNKS();
2854  const int ord0 = kv[0]->GetOrder();
2855  const int nks1 = kv[1]->GetNKS();
2856  const int ord1 = kv[1]->GetOrder();
2857  for (int j = 0; j < nks1; j++)
2858  {
2859  if (kv[1]->isElement(j))
2860  {
2861  for (int i = 0; i < nks0; i++)
2862  {
2863  if (kv[0]->isElement(i))
2864  {
2865  if (activeBdrElem[gbe])
2866  {
2867  Connection conn(lbe,0);
2868  for (int jj = 0; jj <= ord1; jj++)
2869  {
2870  const int jj_ = (okv[1] >= 0) ? (j+jj) : (ny-j-jj);
2871  for (int ii = 0; ii <= ord0; ii++)
2872  {
2873  const int ii_ = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
2874  conn.to = DofMap(p2g(ii_, jj_));
2875  bel_dof_list.Append(conn);
2876  }
2877  }
2878  bel_to_patch[lbe] = b;
2879  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2880  bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
2881  lbe++;
2882  }
2883  gbe++;
2884  }
2885  }
2886  }
2887  }
2888  }
2889  // We must NOT sort bel_dof_list in this case.
2890  bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
2891 }
2892 
2894 {
2895  lvert_vert.SetSize(GetNV());
2896  for (int gv = 0; gv < GetGNV(); gv++)
2897  if (activeVert[gv] >= 0)
2898  {
2899  lvert_vert[activeVert[gv]] = gv;
2900  }
2901 }
2902 
2904 {
2905  lelem_elem.SetSize(GetNE());
2906  for (int le = 0, ge = 0; ge < GetGNE(); ge++)
2907  if (activeElem[ge])
2908  {
2909  lelem_elem[le++] = ge;
2910  }
2911 }
2912 
2913 void NURBSExtension::LoadFE(int i, const FiniteElement *FE) const
2914 {
2915  const NURBSFiniteElement *NURBSFE =
2916  dynamic_cast<const NURBSFiniteElement *>(FE);
2917 
2918  if (NURBSFE->GetElement() != i)
2919  {
2920  Array<int> dofs;
2921  NURBSFE->SetIJK(el_to_IJK.GetRow(i));
2922  if (el_to_patch[i] != NURBSFE->GetPatch())
2923  {
2924  GetPatchKnotVectors(el_to_patch[i], NURBSFE->KnotVectors());
2925  NURBSFE->SetPatch(el_to_patch[i]);
2926  NURBSFE->SetOrder();
2927  }
2928  el_dof->GetRow(i, dofs);
2929  weights.GetSubVector(dofs, NURBSFE->Weights());
2930  NURBSFE->SetElement(i);
2931  }
2932 }
2933 
2934 void NURBSExtension::LoadBE(int i, const FiniteElement *BE) const
2935 {
2936  const NURBSFiniteElement *NURBSFE =
2937  dynamic_cast<const NURBSFiniteElement *>(BE);
2938 
2939  if (NURBSFE->GetElement() != i)
2940  {
2941  Array<int> dofs;
2942  NURBSFE->SetIJK(bel_to_IJK.GetRow(i));
2943  if (bel_to_patch[i] != NURBSFE->GetPatch())
2944  {
2945  GetBdrPatchKnotVectors(bel_to_patch[i], NURBSFE->KnotVectors());
2946  NURBSFE->SetPatch(bel_to_patch[i]);
2947  NURBSFE->SetOrder();
2948  }
2949  bel_dof->GetRow(i, dofs);
2950  weights.GetSubVector(dofs, NURBSFE->Weights());
2951  NURBSFE->SetElement(i);
2952  }
2953 }
2954 
2956 {
2957  delete el_dof;
2958  delete bel_dof;
2959 
2960  if (patches.Size() == 0)
2961  {
2962  GetPatchNets(Nodes, Dimension());
2963  }
2964 }
2965 
2967 {
2968  if (patches.Size() == 0) { return; }
2969 
2970  SetSolutionVector(Nodes, Dimension());
2971  patches.SetSize(0);
2972 }
2973 
2975 {
2976  if (patches.Size() == 0)
2977  {
2978  mfem_error("NURBSExtension::SetKnotsFromPatches :"
2979  " No patches available!");
2980  }
2981 
2983 
2984  for (int p = 0; p < patches.Size(); p++)
2985  {
2986  GetPatchKnotVectors(p, kv);
2987 
2988  for (int i = 0; i < kv.Size(); i++)
2989  {
2990  *kv[i] = *patches[p]->GetKV(i);
2991  }
2992  }
2993 
2994  SetOrdersFromKnotVectors();
2995 
2996  GenerateOffsets();
2997  CountElements();
2998  CountBdrElements();
2999 
3000  // all elements must be active
3001  NumOfActiveElems = NumOfElements;
3002  activeElem.SetSize(NumOfElements);
3003  activeElem = true;
3004 
3005  GenerateActiveVertices();
3006  InitDofMap();
3007  GenerateElementDofTable();
3008  GenerateActiveBdrElems();
3009  GenerateBdrElementDofTable();
3010 
3011  ConnectBoundaries();
3012 }
3013 
3014 void NURBSExtension::LoadSolution(std::istream &input, GridFunction &sol) const
3015 {
3016  const FiniteElementSpace *fes = sol.FESpace();
3017  MFEM_VERIFY(fes->GetNURBSext() == this, "");
3018 
3019  sol.SetSize(fes->GetVSize());
3020 
3021  Array<const KnotVector *> kv(Dimension());
3022  NURBSPatchMap p2g(this);
3023  const int vdim = fes->GetVDim();
3024 
3025  for (int p = 0; p < GetNP(); p++)
3026  {
3027  skip_comment_lines(input, '#');
3028 
3029  p2g.SetPatchDofMap(p, kv);
3030  const int nx = kv[0]->GetNCP();
3031  const int ny = kv[1]->GetNCP();
3032  const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
3033  for (int k = 0; k < nz; k++)
3034  {
3035  for (int j = 0; j < ny; j++)
3036  {
3037  for (int i = 0; i < nx; i++)
3038  {
3039  const int ll = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
3040  const int l = DofMap(ll);
3041  for (int vd = 0; vd < vdim; vd++)
3042  {
3043  input >> sol(fes->DofToVDof(l,vd));
3044  }
3045  }
3046  }
3047  }
3048  }
3049 }
3050 
3051 void NURBSExtension::PrintSolution(const GridFunction &sol, std::ostream &os)
3052 const
3053 {
3054  const FiniteElementSpace *fes = sol.FESpace();
3055  MFEM_VERIFY(fes->GetNURBSext() == this, "");
3056 
3057  Array<const KnotVector *> kv(Dimension());
3058  NURBSPatchMap p2g(this);
3059  const int vdim = fes->GetVDim();
3060 
3061  for (int p = 0; p < GetNP(); p++)
3062  {
3063  os << "\n# patch " << p << "\n\n";
3064 
3065  p2g.SetPatchDofMap(p, kv);
3066  const int nx = kv[0]->GetNCP();
3067  const int ny = kv[1]->GetNCP();
3068  const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
3069  for (int k = 0; k < nz; k++)
3070  {
3071  for (int j = 0; j < ny; j++)
3072  {
3073  for (int i = 0; i < nx; i++)
3074  {
3075  const int ll = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
3076  const int l = DofMap(ll);
3077  os << sol(fes->DofToVDof(l,0));
3078  for (int vd = 1; vd < vdim; vd++)
3079  {
3080  os << ' ' << sol(fes->DofToVDof(l,vd));
3081  }
3082  os << '\n';
3083  }
3084  }
3085  }
3086  }
3087 }
3088 
3089 void NURBSExtension::DegreeElevate(int rel_degree, int degree)
3090 {
3091  for (int p = 0; p < patches.Size(); p++)
3092  {
3093  for (int dir = 0; dir < patches[p]->GetNKV(); dir++)
3094  {
3095  int oldd = patches[p]->GetKV(dir)->GetOrder();
3096  int newd = std::min(oldd + rel_degree, degree);
3097  if (newd > oldd)
3098  {
3099  patches[p]->DegreeElevate(dir, newd - oldd);
3100  }
3101  }
3102  }
3103 }
3104 
3106 {
3107  for (int p = 0; p < patches.Size(); p++)
3108  {
3109  patches[p]->UniformRefinement();
3110  }
3111 }
3112 
3114 {
3115  Array<int> edges;
3116  Array<int> orient;
3117 
3118  Array<KnotVector *> pkv(Dimension());
3119 
3120  for (int p = 0; p < patches.Size(); p++)
3121  {
3122  patchTopo->GetElementEdges(p, edges, orient);
3123 
3124  if (Dimension()==2)
3125  {
3126  pkv[0] = kv[KnotInd(edges[0])];
3127  pkv[1] = kv[KnotInd(edges[1])];
3128  }
3129  else
3130  {
3131  pkv[0] = kv[KnotInd(edges[0])];
3132  pkv[1] = kv[KnotInd(edges[3])];
3133  pkv[2] = kv[KnotInd(edges[8])];
3134  }
3135 
3136  patches[p]->KnotInsert(pkv);
3137  }
3138 }
3139 
3141 {
3142  Array<int> edges;
3143  Array<int> orient;
3144 
3145  Array<Vector *> pkv(Dimension());
3146 
3147  for (int p = 0; p < patches.Size(); p++)
3148  {
3149  patchTopo->GetElementEdges(p, edges, orient);
3150 
3151  if (Dimension()==2)
3152  {
3153  pkv[0] = kv[KnotInd(edges[0])];
3154  pkv[1] = kv[KnotInd(edges[1])];
3155  }
3156  else
3157  {
3158  pkv[0] = kv[KnotInd(edges[0])];
3159  pkv[1] = kv[KnotInd(edges[3])];
3160  pkv[2] = kv[KnotInd(edges[8])];
3161  }
3162 
3163  patches[p]->KnotInsert(pkv);
3164  }
3165 }
3166 
3167 
3168 void NURBSExtension::GetPatchNets(const Vector &coords, int vdim)
3169 {
3170  if (Dimension() == 2)
3171  {
3172  Get2DPatchNets(coords, vdim);
3173  }
3174  else
3175  {
3176  Get3DPatchNets(coords, vdim);
3177  }
3178 }
3179 
3180 void NURBSExtension::Get2DPatchNets(const Vector &coords, int vdim)
3181 {
3183  NURBSPatchMap p2g(this);
3184 
3185  patches.SetSize(GetNP());
3186  for (int p = 0; p < GetNP(); p++)
3187  {
3188  p2g.SetPatchDofMap(p, kv);
3189  patches[p] = new NURBSPatch(kv, vdim+1);
3190  NURBSPatch &Patch = *patches[p];
3191 
3192  for (int j = 0; j < kv[1]->GetNCP(); j++)
3193  {
3194  for (int i = 0; i < kv[0]->GetNCP(); i++)
3195  {
3196  const int l = DofMap(p2g(i,j));
3197  for (int d = 0; d < vdim; d++)
3198  {
3199  Patch(i,j,d) = coords(l*vdim + d)*weights(l);
3200  }
3201  Patch(i,j,vdim) = weights(l);
3202  }
3203  }
3204  }
3205 }
3206 
3207 void NURBSExtension::Get3DPatchNets(const Vector &coords, int vdim)
3208 {
3210  NURBSPatchMap p2g(this);
3211 
3212  patches.SetSize(GetNP());
3213  for (int p = 0; p < GetNP(); p++)
3214  {
3215  p2g.SetPatchDofMap(p, kv);
3216  patches[p] = new NURBSPatch(kv, vdim+1);
3217  NURBSPatch &Patch = *patches[p];
3218 
3219  for (int k = 0; k < kv[2]->GetNCP(); k++)
3220  {
3221  for (int j = 0; j < kv[1]->GetNCP(); j++)
3222  {
3223  for (int i = 0; i < kv[0]->GetNCP(); i++)
3224  {
3225  const int l = DofMap(p2g(i,j,k));
3226  for (int d = 0; d < vdim; d++)
3227  {
3228  Patch(i,j,k,d) = coords(l*vdim + d)*weights(l);
3229  }
3230  Patch(i,j,k,vdim) = weights(l);
3231  }
3232  }
3233  }
3234  }
3235 }
3236 
3238 {
3239  if (Dimension() == 2)
3240  {
3241  Set2DSolutionVector(coords, vdim);
3242  }
3243  else
3244  {
3245  Set3DSolutionVector(coords, vdim);
3246  }
3247 }
3248 
3250 {
3252  NURBSPatchMap p2g(this);
3253 
3254  weights.SetSize(GetNDof());
3255  for (int p = 0; p < GetNP(); p++)
3256  {
3257  p2g.SetPatchDofMap(p, kv);
3258  NURBSPatch &Patch = *patches[p];
3259  MFEM_ASSERT(vdim+1 == Patch.GetNC(), "");
3260 
3261  for (int j = 0; j < kv[1]->GetNCP(); j++)
3262  {
3263  for (int i = 0; i < kv[0]->GetNCP(); i++)
3264  {
3265  const int l = p2g(i,j);
3266  for (int d = 0; d < vdim; d++)
3267  {
3268  coords(l*vdim + d) = Patch(i,j,d)/Patch(i,j,vdim);
3269  }
3270  weights(l) = Patch(i,j,vdim);
3271  }
3272  }
3273  delete patches[p];
3274  }
3275 }
3276 
3278 {
3280  NURBSPatchMap p2g(this);
3281 
3282  weights.SetSize(GetNDof());
3283  for (int p = 0; p < GetNP(); p++)
3284  {
3285  p2g.SetPatchDofMap(p, kv);
3286  NURBSPatch &Patch = *patches[p];
3287  MFEM_ASSERT(vdim+1 == Patch.GetNC(), "");
3288 
3289  for (int k = 0; k < kv[2]->GetNCP(); k++)
3290  {
3291  for (int j = 0; j < kv[1]->GetNCP(); j++)
3292  {
3293  for (int i = 0; i < kv[0]->GetNCP(); i++)
3294  {
3295  const int l = p2g(i,j,k);
3296  for (int d = 0; d < vdim; d++)
3297  {
3298  coords(l*vdim + d) = Patch(i,j,k,d)/Patch(i,j,k,vdim);
3299  }
3300  weights(l) = Patch(i,j,k,vdim);
3301  }
3302  }
3303  }
3304  delete patches[p];
3305  }
3306 }
3307 
3308 
3309 #ifdef MFEM_USE_MPI
3311  : NURBSExtension(orig),
3312  partitioning(orig.partitioning ? new int[orig.GetGNE()] : NULL),
3313  gtopo(orig.gtopo),
3314  ldof_group(orig.ldof_group)
3315 {
3316  // Copy the partitioning, if not NULL
3317  if (partitioning)
3318  {
3319  std::memcpy(partitioning, orig.partitioning, orig.GetGNE()*sizeof(int));
3320  }
3321 }
3322 
3324  int *part, const Array<bool> &active_bel)
3325  : gtopo(comm)
3326 {
3327  if (parent->NumOfActiveElems < parent->NumOfElements)
3328  {
3329  // SetActive (BuildGroups?) and the way the weights are copied
3330  // do not support this case
3331  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
3332  " all elements in the parent must be active!");
3333  }
3334 
3335  patchTopo = parent->patchTopo;
3336  // steal ownership of patchTopo from the 'parent' NURBS extension
3337  if (!parent->own_topo)
3338  {
3339  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
3340  " parent does not own the patch topology!");
3341  }
3342  own_topo = 1;
3343  parent->own_topo = 0;
3344 
3345  parent->edge_to_knot.Copy(edge_to_knot);
3346 
3347  parent->GetOrders().Copy(mOrders);
3348  mOrder = parent->GetOrder();
3349 
3350  NumOfKnotVectors = parent->GetNKV();
3351  knotVectors.SetSize(NumOfKnotVectors);
3352  for (int i = 0; i < NumOfKnotVectors; i++)
3353  {
3354  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
3355  }
3356 
3357  GenerateOffsets();
3358  CountElements();
3359  CountBdrElements();
3360 
3361  // copy 'part' to 'partitioning'
3362  partitioning = new int[GetGNE()];
3363  for (int i = 0; i < GetGNE(); i++)
3364  {
3365  partitioning[i] = part[i];
3366  }
3367  SetActive(partitioning, active_bel);
3368 
3371  // GenerateActiveBdrElems(); // done by SetActive for now
3373 
3374  Table *serial_elem_dof = parent->GetElementDofTable();
3375  BuildGroups(partitioning, *serial_elem_dof);
3376 
3377  weights.SetSize(GetNDof());
3378  // copy weights from parent
3379  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
3380  {
3381  if (activeElem[gel])
3382  {
3383  int ndofs = el_dof->RowSize(lel);
3384  int *ldofs = el_dof->GetRow(lel);
3385  int *gdofs = serial_elem_dof->GetRow(gel);
3386  for (int i = 0; i < ndofs; i++)
3387  {
3388  weights(ldofs[i]) = parent->weights(gdofs[i]);
3389  }
3390  lel++;
3391  }
3392  }
3393 }
3394 
3396  const ParNURBSExtension *par_parent)
3397  : gtopo(par_parent->gtopo.GetComm())
3398 {
3399  // steal all data from parent
3400  mOrder = parent->mOrder;
3401  Swap(mOrders, parent->mOrders);
3402 
3403  patchTopo = parent->patchTopo;
3404  own_topo = parent->own_topo;
3405  parent->own_topo = 0;
3406 
3407  Swap(edge_to_knot, parent->edge_to_knot);
3408 
3410  Swap(knotVectors, parent->knotVectors);
3411 
3412  NumOfVertices = parent->NumOfVertices;
3413  NumOfElements = parent->NumOfElements;
3415  NumOfDofs = parent->NumOfDofs;
3416 
3417  Swap(v_meshOffsets, parent->v_meshOffsets);
3418  Swap(e_meshOffsets, parent->e_meshOffsets);
3419  Swap(f_meshOffsets, parent->f_meshOffsets);
3420  Swap(p_meshOffsets, parent->p_meshOffsets);
3421 
3426 
3427  Swap(d_to_d, parent->d_to_d);
3428  Swap(master, parent->master);
3429  Swap(slave, parent->slave);
3430 
3434  NumOfActiveDofs = parent->NumOfActiveDofs;
3435 
3436  Swap(activeVert, parent->activeVert);
3437  Swap(activeElem, parent->activeElem);
3438  Swap(activeBdrElem, parent->activeBdrElem);
3439  Swap(activeDof, parent->activeDof);
3440 
3441  el_dof = parent->el_dof;
3442  bel_dof = parent->bel_dof;
3443  parent->el_dof = parent->bel_dof = NULL;
3444 
3445  Swap(el_to_patch, parent->el_to_patch);
3446  Swap(bel_to_patch, parent->bel_to_patch);
3447  Swap(el_to_IJK, parent->el_to_IJK);
3448  Swap(bel_to_IJK, parent->bel_to_IJK);
3449 
3450  Swap(weights, parent->weights);
3451  MFEM_VERIFY(!parent->HavePatches(), "");
3452 
3453  delete parent;
3454 
3455  partitioning = NULL;
3456 
3457  MFEM_VERIFY(par_parent->partitioning,
3458  "parent ParNURBSExtension has no partitioning!");
3459 
3460  // Support for the case when 'parent' is not a local NURBSExtension, i.e.
3461  // NumOfActiveElems is not the same as in 'par_parent'. In that case, we
3462  // assume 'parent' is a global NURBSExtension, i.e. all elements are active.
3463  bool extract_weights = false;
3464  if (NumOfActiveElems != par_parent->NumOfActiveElems)
3465  {
3466  MFEM_ASSERT(NumOfActiveElems == NumOfElements, "internal error");
3467 
3468  SetActive(par_parent->partitioning, par_parent->activeBdrElem);
3470  delete el_dof;
3472  el_to_IJK.DeleteAll();
3474  // GenerateActiveBdrElems(); // done by SetActive for now
3475  delete bel_dof;
3479  extract_weights = true;
3480  }
3481 
3482  Table *glob_elem_dof = GetGlobalElementDofTable();
3483  BuildGroups(par_parent->partitioning, *glob_elem_dof);
3484  if (extract_weights)
3485  {
3486  Vector glob_weights;
3487  Swap(weights, glob_weights);
3488  weights.SetSize(GetNDof());
3489  // Copy the local 'weights' from the 'glob_weights'.
3490  // Assumption: the local element ids follow the global ordering.
3491  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
3492  {
3493  if (activeElem[gel])
3494  {
3495  int ndofs = el_dof->RowSize(lel);
3496  int *ldofs = el_dof->GetRow(lel);
3497  int *gdofs = glob_elem_dof->GetRow(gel);
3498  for (int i = 0; i < ndofs; i++)
3499  {
3500  weights(ldofs[i]) = glob_weights(gdofs[i]);
3501  }
3502  lel++;
3503  }
3504  }
3505  }
3506  delete glob_elem_dof;
3507 }
3508 
3509 Table *ParNURBSExtension::GetGlobalElementDofTable()
3510 {
3511  if (Dimension() == 2)
3512  {
3513  return Get2DGlobalElementDofTable();
3514  }
3515  else
3516  {
3517  return Get3DGlobalElementDofTable();
3518  }
3519 }
3520 
3521 Table *ParNURBSExtension::Get2DGlobalElementDofTable()
3522 {
3523  int el = 0;
3524  const KnotVector *kv[2];
3525  NURBSPatchMap p2g(this);
3526  Array<Connection> gel_dof_list;
3527 
3528  for (int p = 0; p < GetNP(); p++)
3529  {
3530  p2g.SetPatchDofMap(p, kv);
3531 
3532  // Load dofs
3533  const int ord0 = kv[0]->GetOrder();
3534  const int ord1 = kv[1]->GetOrder();
3535  for (int j = 0; j < kv[1]->GetNKS(); j++)
3536  {
3537  if (kv[1]->isElement(j))
3538  {
3539  for (int i = 0; i < kv[0]->GetNKS(); i++)
3540  {
3541  if (kv[0]->isElement(i))
3542  {
3543  Connection conn(el,0);
3544  for (int jj = 0; jj <= ord1; jj++)
3545  {
3546  for (int ii = 0; ii <= ord0; ii++)
3547  {
3548  conn.to = p2g(i+ii,j+jj);
3549  gel_dof_list.Append(conn);
3550  }
3551  }
3552  el++;
3553  }
3554  }
3555  }
3556  }
3557  }
3558  // We must NOT sort gel_dof_list in this case.
3559  return (new Table(GetGNE(), gel_dof_list));
3560 }
3561 
3562 Table *ParNURBSExtension::Get3DGlobalElementDofTable()
3563 {
3564  int el = 0;
3565  const KnotVector *kv[3];
3566  NURBSPatchMap p2g(this);
3567  Array<Connection> gel_dof_list;
3568 
3569  for (int p = 0; p < GetNP(); p++)
3570  {
3571  p2g.SetPatchDofMap(p, kv);
3572 
3573  // Load dofs
3574  const int ord0 = kv[0]->GetOrder();
3575  const int ord1 = kv[1]->GetOrder();
3576  const int ord2 = kv[2]->GetOrder();
3577  for (int k = 0; k < kv[2]->GetNKS(); k++)
3578  {
3579  if (kv[2]->isElement(k))
3580  {
3581  for (int j = 0; j < kv[1]->GetNKS(); j++)
3582  {
3583  if (kv[1]->isElement(j))
3584  {
3585  for (int i = 0; i < kv[0]->GetNKS(); i++)
3586  {
3587  if (kv[0]->isElement(i))
3588  {
3589  Connection conn(el,0);
3590  for (int kk = 0; kk <= ord2; kk++)
3591  {
3592  for (int jj = 0; jj <= ord1; jj++)
3593  {
3594  for (int ii = 0; ii <= ord0; ii++)
3595  {
3596  conn.to = p2g(i+ii,j+jj,k+kk);
3597  gel_dof_list.Append(conn);
3598  }
3599  }
3600  }
3601  el++;
3602  }
3603  }
3604  }
3605  }
3606  }
3607  }
3608  }
3609  // We must NOT sort gel_dof_list in this case.
3610  return (new Table(GetGNE(), gel_dof_list));
3611 }
3612 
3613 void ParNURBSExtension::SetActive(const int *partitioning_,
3614  const Array<bool> &active_bel)
3615 {
3617  activeElem = false;
3618  NumOfActiveElems = 0;
3619  const int MyRank = gtopo.MyRank();
3620  for (int i = 0; i < GetGNE(); i++)
3621  if (partitioning_[i] == MyRank)
3622  {
3623  activeElem[i] = true;
3624  NumOfActiveElems++;
3625  }
3626 
3627  active_bel.Copy(activeBdrElem);
3628  NumOfActiveBdrElems = 0;
3629  for (int i = 0; i < GetGNBE(); i++)
3630  if (activeBdrElem[i])
3631  {
3633  }
3634 }
3635 
3636 void ParNURBSExtension::BuildGroups(const int *partitioning_,
3637  const Table &elem_dof)
3638 {
3639  Table dof_proc;
3640 
3641  ListOfIntegerSets groups;
3642  IntegerSet group;
3643 
3644  Transpose(elem_dof, dof_proc); // dof_proc is dof_elem
3645  // convert elements to processors
3646  for (int i = 0; i < dof_proc.Size_of_connections(); i++)
3647  {
3648  dof_proc.GetJ()[i] = partitioning_[dof_proc.GetJ()[i]];
3649  }
3650 
3651  // the first group is the local one
3652  int MyRank = gtopo.MyRank();
3653  group.Recreate(1, &MyRank);
3654  groups.Insert(group);
3655 
3656  int dof = 0;
3658  for (int d = 0; d < GetNTotalDof(); d++)
3659  if (activeDof[d])
3660  {
3661  group.Recreate(dof_proc.RowSize(d), dof_proc.GetRow(d));
3662  ldof_group[dof] = groups.Insert(group);
3663 
3664  dof++;
3665  }
3666 
3667  gtopo.Create(groups, 1822);
3668 }
3669 #endif // MFEM_USE_MPI
3670 
3671 
3672 void NURBSPatchMap::GetPatchKnotVectors(int p, const KnotVector *kv[])
3673 {
3674  Ext->patchTopo->GetElementVertices(p, verts);
3675  Ext->patchTopo->GetElementEdges(p, edges, oedge);
3676  if (Ext->Dimension() == 2)
3677  {
3678  kv[0] = Ext->KnotVec(edges[0]);
3679  kv[1] = Ext->KnotVec(edges[1]);
3680  }
3681  else
3682  {
3683  Ext->patchTopo->GetElementFaces(p, faces, oface);
3684 
3685  kv[0] = Ext->KnotVec(edges[0]);
3686  kv[1] = Ext->KnotVec(edges[3]);
3687  kv[2] = Ext->KnotVec(edges[8]);
3688  }
3689  opatch = 0;
3690 }
3691 
3692 void NURBSPatchMap::GetBdrPatchKnotVectors(int p, const KnotVector *kv[],
3693  int *okv)
3694 {
3695  Ext->patchTopo->GetBdrElementVertices(p, verts);
3696  Ext->patchTopo->GetBdrElementEdges(p, edges, oedge);
3697  kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]);
3698  if (Ext->Dimension() == 3)
3699  {
3700  faces.SetSize(1);
3701  Ext->patchTopo->GetBdrElementFace(p, &faces[0], &opatch);
3702  kv[1] = Ext->KnotVec(edges[1], oedge[1], &okv[1]);
3703  }
3704  else
3705  {
3706  opatch = oedge[0];
3707  }
3708 }
3709 
3711 {
3712  GetPatchKnotVectors(p, kv);
3713 
3714  I = kv[0]->GetNE() - 1;
3715  J = kv[1]->GetNE() - 1;
3716 
3717  for (int i = 0; i < verts.Size(); i++)
3718  {
3719  verts[i] = Ext->v_meshOffsets[verts[i]];
3720  }
3721 
3722  for (int i = 0; i < edges.Size(); i++)
3723  {
3724  edges[i] = Ext->e_meshOffsets[edges[i]];
3725  }
3726 
3727  if (Ext->Dimension() == 3)
3728  {
3729  K = kv[2]->GetNE() - 1;
3730 
3731  for (int i = 0; i < faces.Size(); i++)
3732  {
3733  faces[i] = Ext->f_meshOffsets[faces[i]];
3734  }
3735  }
3736 
3737  pOffset = Ext->p_meshOffsets[p];
3738 }
3739 
3741 {
3742  GetPatchKnotVectors(p, kv);
3743 
3744  I = kv[0]->GetNCP() - 2;
3745  J = kv[1]->GetNCP() - 2;
3746 
3747  for (int i = 0; i < verts.Size(); i++)
3748  {
3749  verts[i] = Ext->v_spaceOffsets[verts[i]];
3750  }
3751 
3752  for (int i = 0; i < edges.Size(); i++)
3753  {
3754  edges[i] = Ext->e_spaceOffsets[edges[i]];
3755  }
3756 
3757  if (Ext->Dimension() == 3)
3758  {
3759  K = kv[2]->GetNCP() - 2;
3760 
3761  for (int i = 0; i < faces.Size(); i++)
3762  {
3763  faces[i] = Ext->f_spaceOffsets[faces[i]];
3764  }
3765  }
3766 
3767  pOffset = Ext->p_spaceOffsets[p];
3768 }
3769 
3771  int *okv)
3772 {
3773  GetBdrPatchKnotVectors(p, kv, okv);
3774 
3775  I = kv[0]->GetNE() - 1;
3776 
3777  for (int i = 0; i < verts.Size(); i++)
3778  {
3779  verts[i] = Ext->v_meshOffsets[verts[i]];
3780  }
3781 
3782  if (Ext->Dimension() == 2)
3783  {
3784  pOffset = Ext->e_meshOffsets[edges[0]];
3785  }
3786  else
3787  {
3788  J = kv[1]->GetNE() - 1;
3789 
3790  for (int i = 0; i < edges.Size(); i++)
3791  {
3792  edges[i] = Ext->e_meshOffsets[edges[i]];
3793  }
3794 
3795  pOffset = Ext->f_meshOffsets[faces[0]];
3796  }
3797 }
3798 
3799 void NURBSPatchMap::SetBdrPatchDofMap(int p, const KnotVector *kv[], int *okv)
3800 {
3801  GetBdrPatchKnotVectors(p, kv, okv);
3802 
3803  I = kv[0]->GetNCP() - 2;
3804 
3805  for (int i = 0; i < verts.Size(); i++)
3806  {
3807  verts[i] = Ext->v_spaceOffsets[verts[i]];
3808  }
3809 
3810  if (Ext->Dimension() == 2)
3811  {
3812  pOffset = Ext->e_spaceOffsets[edges[0]];
3813  }
3814  else
3815  {
3816  J = kv[1]->GetNCP() - 2;
3817 
3818  for (int i = 0; i < edges.Size(); i++)
3819  {
3820  edges[i] = Ext->e_spaceOffsets[edges[i]];
3821  }
3822 
3823  pOffset = Ext->f_spaceOffsets[faces[0]];
3824  }
3825 }
3826 
3827 }
Array< KnotVector * > knotVectors
Definition: nurbs.hpp:191
Abstract class for all finite elements.
Definition: fe_base.hpp:235
Array< int > f_meshOffsets
Definition: nurbs.hpp:204
KnotVector * DegreeElevate(int t) const
Definition: nurbs.cpp:58
void Create(ListOfIntegerSets &groups, int mpitag)
Set up the group topology given the list of sets of shared entities.
Array2D< int > bel_to_IJK
Definition: nurbs.hpp:218
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:573
void SetCoordsFromPatches(Vector &Nodes)
Definition: nurbs.cpp:2966
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
int GetGNE() const
Definition: nurbs.hpp:359
Array< int > activeVert
Definition: nurbs.hpp:183
void ConnectBoundaries2D(int bnd0, int bnd1)
Definition: nurbs.cpp:1824
KnotVector * KnotVec(int edge)
Definition: nurbs.hpp:557
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:251
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
Definition: mesh.cpp:5954
void DegreeElevate(int dir, int t)
Definition: nurbs.cpp:837
int NumOfElements
Definition: nurbs.hpp:37
NURBSPatch * Interpolate(NURBSPatch &p1, NURBSPatch &p2)
Definition: nurbs.cpp:1213
void LoadBE(int i, const FiniteElement *BE) const
Definition: nurbs.cpp:2934
void Generate2DElementDofTable()
Definition: nurbs.cpp:2653
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Definition: nurbs.cpp:2115
Array< int > ldof_group
Definition: nurbs.hpp:425
int GetPatch() const
Definition: fe_nurbs.hpp:51
Array< int > master
Definition: nurbs.hpp:198
Array< const KnotVector * > & KnotVectors() const
Definition: fe_nurbs.hpp:55
void Generate3DElementDofTable()
Definition: nurbs.cpp:2707
void GetBdrElementFace(int i, int *f, int *o) const
Return the index and the orientation of the face of bdr element i. (3D)
Definition: mesh.cpp:6171
void KnotInsert(int dir, const KnotVector &knot)
Definition: nurbs.cpp:689
void MergeWeights(Mesh *mesh_array[], int num_pieces)
Definition: nurbs.cpp:2090
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
void Get3DBdrElementTopo(Array< Element * > &boundary) const
Definition: nurbs.cpp:2585
void UniformRefinement()
Definition: nurbs.cpp:3105
int GetNE() const
Definition: nurbs.hpp:48
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:856
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:547
void Rotate3D(double normal[], double angle)
Definition: nurbs.cpp:1164
Array< int > e_meshOffsets
Definition: nurbs.hpp:203
int GetNKV() const
Definition: nurbs.hpp:355
int GetNKS() const
Definition: nurbs.hpp:49
void GetElements()
Count the number of elements.
Definition: nurbs.cpp:101
void SetOrdersFromKnotVectors()
Definition: nurbs.cpp:2308
An arbitrary order and dimension NURBS element.
Definition: fe_nurbs.hpp:23
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
void Generate3DBdrElementDofTable()
Definition: nurbs.cpp:2835
int NumOfControlPoints
Definition: nurbs.hpp:37
Array< int > e_spaceOffsets
Definition: nurbs.hpp:209
void GenerateBdrElementDofTable()
Definition: nurbs.cpp:2773
const KnotVector * GetKnotVector(int i) const
Definition: nurbs.hpp:368
void GetElementLocalToGlobal(Array< int > &lelem_elem)
Definition: nurbs.cpp:2903
void init(int dim_)
Definition: nurbs.cpp:407
Data type quadrilateral element.
Array< int > edge_to_knot
Definition: nurbs.hpp:190
void GenerateActiveVertices()
Definition: nurbs.cpp:2001
Array< bool > activeElem
Definition: nurbs.hpp:184
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
Definition: nurbs.cpp:2893
void Set2DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:3249
Vector knot
Definition: nurbs.hpp:36
int GetElement() const
Definition: fe_nurbs.hpp:53
int GetNTotalDof() const
Definition: nurbs.hpp:364
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
Definition: text.hpp:31
void DeleteAll()
Delete the whole array.
Definition: array.hpp:846
int GetNDof() const
Definition: nurbs.hpp:365
int Dimension() const
Definition: nurbs.hpp:345
Array< int > el_to_patch
Definition: nurbs.hpp:215
void SwapDirections(int dir1, int dir2)
Definition: nurbs.cpp:1100
friend class NURBSPatchMap
Definition: nurbs.hpp:171
KnotVector()
Create KnotVector.
Definition: nurbs.hpp:41
static const int MaxOrder
Definition: nurbs.hpp:34
Array< int > v_spaceOffsets
Definition: nurbs.hpp:208
void UniformRefinement(Vector &newknots) const
Definition: nurbs.cpp:87
void Generate2DBdrElementDofTable()
Definition: nurbs.cpp:2792
Data type hexahedron element.
Definition: hexahedron.hpp:22
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
void ConnectBoundaries3D(int bnd0, int bnd1)
Definition: nurbs.cpp:1900
ParNURBSExtension(const ParNURBSExtension &orig)
Definition: nurbs.cpp:3310
void SetIJK(const int *IJK) const
Definition: fe_nurbs.hpp:50
Array< int > d_to_d
Definition: nurbs.hpp:197
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
void Set3DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:3277
void SetPatch(int p) const
Definition: fe_nurbs.hpp:52
void GetPatchKnotVectors(int p, Array< KnotVector * > &kv)
Definition: nurbs.cpp:2214
void GenerateActiveBdrElems()
Definition: nurbs.cpp:2068
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
void PrintSolution(const GridFunction &sol, std::ostream &out) const
Definition: nurbs.cpp:3051
void Get3DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3207
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
int MyRank() const
Return the MPI rank within this object&#39;s communicator.
double b
Definition: lissajous.cpp:42
void Get2DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3180
bool HavePatches() const
Definition: nurbs.hpp:374
void GenerateOffsets()
Definition: nurbs.cpp:2318
Array< int > activeDof
Definition: nurbs.hpp:186
void GetBdrElementTopo(Array< Element * > &boundary) const
Definition: nurbs.cpp:2540
void SetKnotsFromPatches()
Definition: nurbs.cpp:2974
void GetElementTopo(Array< Element * > &elements) const
Definition: nurbs.cpp:2445
const Array< int > & GetOrders() const
Read-only access to the orders of all knot vectors.
Definition: nurbs.hpp:350
Array< int > f_spaceOffsets
Definition: nurbs.hpp:210
Array< bool > activeBdrElem
Definition: nurbs.hpp:185
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
int GetNC() const
Definition: nurbs.hpp:136
void LoadSolution(std::istream &input, GridFunction &sol) const
Definition: nurbs.cpp:3014
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:668
void SetElement(int e) const
Definition: fe_nurbs.hpp:54
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:353
Array< int > slave
Definition: nurbs.hpp:199
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:630
void LoadFE(int i, const FiniteElement *FE) const
Definition: nurbs.cpp:2913
const NURBSExtension * GetNURBSext() const
Definition: fespace.hpp:443
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
Array< int > v_meshOffsets
Definition: nurbs.hpp:202
int GetNP() const
Definition: nurbs.hpp:346
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
void SetBdrPatchVertexMap(int p, const KnotVector *kv[], int *okv)
Definition: nurbs.cpp:3770
void GetBdrPatchKnotVectors(int p, Array< KnotVector * > &kv)
Definition: nurbs.cpp:2255
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
Array< int > p_meshOffsets
Definition: nurbs.hpp:205
KnotVector & operator=(const KnotVector &kv)
Definition: nurbs.cpp:46
void SetPatchVertexMap(int p, const KnotVector *kv[])
Definition: nurbs.cpp:3710
Helper struct for defining a connectivity table, see Table::MakeFromList.
Definition: table.hpp:27
Table * GetElementDofTable()
Definition: nurbs.hpp:376
Array2D< int > el_to_IJK
Definition: nurbs.hpp:217
void SetOrderFromOrders()
Definition: nurbs.cpp:2294
string direction
void CalcShape(Vector &shape, int i, double xi) const
Definition: nurbs.cpp:160
void CalcDnShape(Vector &gradn, int n, int i, double xi) const
Definition: nurbs.cpp:243
void CountBdrElements()
Definition: nurbs.cpp:2425
void PrintFunctions(std::ostream &out, int samples=11) const
Definition: nurbs.cpp:133
int GetNE() const
Definition: nurbs.hpp:360
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
void GetPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3168
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
Definition: mesh.cpp:5932
NURBSPatch * Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
Definition: nurbs.cpp:1260
Array< int > p_spaceOffsets
Definition: nurbs.hpp:211
void SetSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:3237
void Get3DElementTopo(Array< Element * > &elements) const
Definition: nurbs.cpp:2495
void SetBdrPatchDofMap(int p, const KnotVector *kv[], int *okv)
Definition: nurbs.cpp:3799
void ConnectBoundaries()
Definition: nurbs.cpp:1793
int GetNCP() const
Definition: nurbs.hpp:50
int dim
Definition: ex24.cpp:53
void GenerateElementDofTable()
Definition: nurbs.cpp:2623
void Difference(const KnotVector &kv, Vector &diff) const
Definition: nurbs.cpp:374
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition: hooke.cpp:45
Vector & Weights() const
Definition: fe_nurbs.hpp:56
Array< int > mOrders
Definition: nurbs.hpp:175
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition: fe_nurbs.hpp:58
RefCoord t[3]
Vector data type.
Definition: vector.hpp:60
void Print(std::ostream &out) const
Definition: nurbs.cpp:126
void DeleteAll()
Delete all dynamically allocated memory, resetting all dimensions to zero.
Definition: array.hpp:439
int GetGNBE() const
Definition: nurbs.hpp:361
int RowSize(int i) const
Definition: table.hpp:108
void CalcDShape(Vector &grad, int i, double xi) const
Definition: nurbs.cpp:186
void DegreeElevate(int rel_degree, int degree=16)
Definition: nurbs.cpp:3089
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
void CheckBdrPatches()
Definition: nurbs.cpp:2186
void Get2DElementTopo(Array< Element * > &elements) const
Definition: nurbs.cpp:2459
GroupTopology gtopo
Definition: nurbs.hpp:423
Data type line segment element.
Definition: segment.hpp:22
int GetOrder() const
Definition: nurbs.hpp:51
Array< int > bel_to_patch
Definition: nurbs.hpp:216
void SetPatchDofMap(int p, const KnotVector *kv[])
Definition: nurbs.cpp:3740
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Definition: mesh.hpp:1113
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:1109
void ConvertToPatches(const Vector &Nodes)
Definition: nurbs.cpp:2955
void KnotInsert(Array< KnotVector * > &kv)
Definition: nurbs.cpp:3113
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
Definition: mesh.cpp:6148
void Get2DBdrElementTopo(Array< Element * > &boundary) const
Definition: nurbs.cpp:2554
int SetLoopDirection(int dir)
Definition: nurbs.cpp:612