MFEM  v3.2
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "../fem/fem.hpp"
13 #include <algorithm>
14 #if defined(_MSC_VER) && (_MSC_VER < 1800)
15 #include <float.h>
16 #define copysign _copysign
17 #endif
18 
19 namespace mfem
20 {
21 
22 using namespace std;
23 
24 const int KnotVector::MaxOrder = 10;
25 
26 KnotVector::KnotVector(std::istream &input)
27 {
28  input >> Order >> NumOfControlPoints;
29  knot.Load(input, NumOfControlPoints + Order + 1);
30  GetElements();
31 }
32 
33 KnotVector::KnotVector(int Order_, int NCP)
34 {
35  Order = Order_;
36  NumOfControlPoints = NCP;
37  knot.SetSize(NumOfControlPoints + Order + 1);
38 
39  knot = -1.;
40 }
41 
43 {
44  Order = kv.Order;
45  NumOfControlPoints = kv.NumOfControlPoints;
46  NumOfElements = kv.NumOfElements;
47  knot = kv.knot;
48  // alternatively, re-compute NumOfElements
49  // GetElements();
50 
51  return *this;
52 }
53 
55 {
56  if (t < 0)
57  mfem_error("KnotVector::DegreeElevate :\n"
58  " Parent KnotVector order higher than child");
59 
60  int nOrder = Order + t;
61  KnotVector *newkv = new KnotVector(nOrder, GetNCP() + t);
62 
63  for (int i = 0; i <= nOrder; i++)
64  {
65  (*newkv)[i] = knot(0);
66  }
67  for (int i = nOrder + 1; i < newkv->GetNCP(); i++)
68  {
69  (*newkv)[i] = knot(i - t);
70  }
71  for (int i = 0; i <= nOrder; i++)
72  {
73  (*newkv)[newkv->GetNCP() + i] = knot(knot.Size()-1);
74  }
75 
76  newkv->GetElements();
77 
78  return newkv;
79 }
80 
82 {
83  newknots.SetSize(NumOfElements);
84  int j = 0;
85  for (int i = 0; i < knot.Size()-1; i++)
86  {
87  if (knot(i) != knot(i+1))
88  {
89  newknots(j) = 0.5*(knot(i) + knot(i+1));
90  j++;
91  }
92  }
93 }
94 
96 {
97  NumOfElements = 0;
98  for (int i = Order; i < NumOfControlPoints; i++)
99  {
100  if (knot(i) != knot(i+1))
101  {
102  NumOfElements++;
103  }
104  }
105 }
106 
108 {
109  double apb = knot(0) + knot(knot.Size()-1);
110 
111  int ns = (NumOfControlPoints - Order)/2;
112  for (int i = 1; i <= ns; i++)
113  {
114  double tmp = apb - knot(Order + i);
115  knot(Order + i) = apb - knot(NumOfControlPoints - i);
116  knot(NumOfControlPoints - i) = tmp;
117  }
118 }
119 
120 void KnotVector::Print(std::ostream &out) const
121 {
122  out << Order << ' ' << NumOfControlPoints << ' ';
123  knot.Print(out, knot.Size());
124 }
125 
126 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
127 void KnotVector::CalcShape(Vector &shape, int i, double xi)
128 {
129  int p = Order;
130  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
131  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), saved, tmp;
132  double left[MaxOrder+1], right[MaxOrder+1];
133 
134 #ifdef MFEM_DEBUG
135  if (p > MaxOrder)
136  {
137  mfem_error("KnotVector::CalcShape : Order > MaxOrder!");
138  }
139 #endif
140 
141  shape(0) = 1.;
142  for (int j = 1; j <= p; ++j)
143  {
144  left[j] = u - knot(ip+1-j);
145  right[j] = knot(ip+j) - u;
146  saved = 0.;
147  for (int r = 0; r < j; ++r)
148  {
149  tmp = shape(r)/(right[r+1] + left[j-r]);
150  shape(r) = saved + right[r+1]*tmp;
151  saved = left[j-r]*tmp;
152  }
153  shape(j) = saved;
154  }
155 }
156 
157 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
158 void KnotVector::CalcDShape(Vector &grad, int i, double xi)
159 {
160  int p = Order, rk, pk;
161  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
162  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), temp, saved, d;
163  double ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1], right[MaxOrder+1];
164 
165 #ifdef MFEM_DEBUG
166  if (p > MaxOrder)
167  {
168  mfem_error("KnotVector::CalcDShape : Order > MaxOrder!");
169  }
170 #endif
171 
172  ndu[0][0] = 1.0;
173  for (int j = 1; j <= p; j++)
174  {
175  left[j] = u - knot(ip-j+1);
176  right[j] = knot(ip+j) - u;
177  saved = 0.0;
178  for (int r = 0; r < j; r++)
179  {
180  ndu[j][r] = right[r+1] + left[j-r];
181  temp = ndu[r][j-1]/ndu[j][r];
182  ndu[r][j] = saved + right[r+1]*temp;
183  saved = left[j-r]*temp;
184  }
185  ndu[j][j] = saved;
186  }
187 
188  for (int r = 0; r <= p; ++r)
189  {
190  d = 0.0;
191  rk = r-1;
192  pk = p-1;
193  if (r >= 1)
194  {
195  d = ndu[rk][pk]/ndu[p][rk];
196  }
197  if (r <= pk)
198  {
199  d -= ndu[r][pk]/ndu[p][r];
200  }
201  grad(r) = d;
202  }
203 
204  if (i >= 0)
205  {
206  grad *= p*(knot(ip+1) - knot(ip));
207  }
208  else
209  {
210  grad *= p*(knot(ip) - knot(ip+1));
211  }
212 }
213 
214 int KnotVector::findKnotSpan(double u) const
215 {
216  int low, mid, high;
217 
218  if (u == knot(NumOfControlPoints+Order))
219  {
220  mid = NumOfControlPoints;
221  }
222  else
223  {
224  low = Order;
225  high = NumOfControlPoints + 1;
226  mid = (low + high)/2;
227  while ( (u < knot(mid-1)) || (u > knot(mid)) )
228  {
229  if (u < knot(mid-1))
230  {
231  high = mid;
232  }
233  else
234  {
235  low = mid;
236  }
237  mid = (low + high)/2;
238  }
239  }
240  return mid;
241 }
242 
243 void KnotVector::Difference(const KnotVector &kv, Vector &diff) const
244 {
245  if (Order != kv.GetOrder())
246  {
247  mfem_error("KnotVector::Difference :\n"
248  " Can not compare knot vectors with different orders!");
249  }
250 
251  int s = kv.Size() - Size();
252  if (s < 0)
253  {
254  kv.Difference(*this, diff);
255  return;
256  }
257 
258  diff.SetSize(s);
259 
260  s = 0;
261  int i = 0;
262  for (int j = 0; j < kv.Size(); j++)
263  {
264  if (knot(i) == kv[j])
265  {
266  i++;
267  }
268  else
269  {
270  diff(s) = kv[j];
271  s++;
272  }
273  }
274 }
275 
276 
277 void NURBSPatch::init(int dim_)
278 {
279  Dim = dim_;
280  sd = nd = -1;
281 
282  if (kv.Size() == 2)
283  {
284  ni = kv[0]->GetNCP();
285  nj = kv[1]->GetNCP();
286  nk = -1;
287 
288  data = new double[ni*nj*Dim];
289 
290 #ifdef MFEM_DEBUG
291  for (int i = 0; i < ni*nj*Dim; i++)
292  {
293  data[i] = -999.99;
294  }
295 #endif
296  }
297  else if (kv.Size() == 3)
298  {
299  ni = kv[0]->GetNCP();
300  nj = kv[1]->GetNCP();
301  nk = kv[2]->GetNCP();
302 
303  data = new double[ni*nj*nk*Dim];
304 
305 #ifdef MFEM_DEBUG
306  for (int i = 0; i < ni*nj*nk*Dim; i++)
307  {
308  data[i] = -999.99;
309  }
310 #endif
311  }
312  else
313  {
314  mfem_error("NURBSPatch::init : Wrond dimension of knotvectors!");
315  }
316 }
317 
318 NURBSPatch::NURBSPatch(std::istream &input)
319 {
320  int pdim, dim, size = 1;
321  string ident;
322 
323  input >> ws >> ident >> pdim; // knotvectors
324  kv.SetSize(pdim);
325  for (int i = 0; i < pdim; i++)
326  {
327  kv[i] = new KnotVector(input);
328  size *= kv[i]->GetNCP();
329  }
330 
331  input >> ws >> ident >> dim; // dimension
332  init(dim + 1);
333 
334  input >> ws >> ident; // controlpoints (homogeneous coordinates)
335  if (ident == "controlpoints" || ident == "controlpoints_homogeneous")
336  {
337  for (int j = 0, i = 0; i < size; i++)
338  for (int d = 0; d <= dim; d++, j++)
339  {
340  input >> data[j];
341  }
342  }
343  else // "controlpoints_cartesian" (Cartesian coordinates with weight)
344  {
345  for (int j = 0, i = 0; i < size; i++)
346  {
347  for (int d = 0; d <= dim; d++)
348  {
349  input >> data[j+d];
350  }
351  for (int d = 0; d < dim; d++)
352  {
353  data[j+d] *= data[j+dim];
354  }
355  j += (dim+1);
356  }
357  }
358 }
359 
361 {
362  kv.SetSize(2);
363  kv[0] = new KnotVector(*kv0);
364  kv[1] = new KnotVector(*kv1);
365  init(dim_);
366 }
367 
369  int dim_)
370 {
371  kv.SetSize(3);
372  kv[0] = new KnotVector(*kv0);
373  kv[1] = new KnotVector(*kv1);
374  kv[2] = new KnotVector(*kv2);
375  init(dim_);
376 }
377 
379 {
380  kv.SetSize(kv_.Size());
381  for (int i = 0; i < kv.Size(); i++)
382  {
383  kv[i] = new KnotVector(*kv_[i]);
384  }
385  init(dim_);
386 }
387 
388 NURBSPatch::NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP)
389 {
390  kv.SetSize(parent->kv.Size());
391  for (int i = 0; i < kv.Size(); i++)
392  if (i != dir)
393  {
394  kv[i] = new KnotVector(*parent->kv[i]);
395  }
396  else
397  {
398  kv[i] = new KnotVector(Order, NCP);
399  }
400  init(parent->Dim);
401 }
402 
404 {
405  if (data != NULL)
406  {
407  delete [] data;
408  }
409 
410  for (int i = 0; i < kv.Size(); i++)
411  {
412  if (kv[i]) { delete kv[i]; }
413  }
414 
415  data = np->data;
416  np->kv.Copy(kv);
417 
418  ni = np->ni;
419  nj = np->nj;
420  nk = np->nk;
421  Dim = np->Dim;
422 
423  np->data = NULL;
424  np->kv.SetSize(0);
425 
426  delete np;
427 }
428 
430 {
431  if (data != NULL)
432  {
433  delete [] data;
434  }
435 
436  for (int i = 0; i < kv.Size(); i++)
437  {
438  if (kv[i]) { delete kv[i]; }
439  }
440 }
441 
442 void NURBSPatch::Print(std::ostream &out)
443 {
444  int size = 1;
445 
446  out << "knotvectors\n" << kv.Size() << '\n';
447  for (int i = 0; i < kv.Size(); i++)
448  {
449  kv[i]->Print(out);
450  size *= kv[i]->GetNCP();
451  }
452 
453  out << "\ndimension\n" << Dim - 1
454  << "\n\ncontrolpoints\n";
455  for (int j = 0, i = 0; i < size; i++)
456  {
457  out << data[j++];
458  for (int d = 1; d < Dim; d++)
459  {
460  out << ' ' << data[j++];
461  }
462  out << '\n';
463  }
464 }
465 
467 {
468  if (nk == -1)
469  {
470  if (dir == 0)
471  {
472  sd = Dim;
473  nd = ni;
474 
475  return nj*Dim;
476  }
477  else if (dir == 1)
478  {
479  sd = ni*Dim;
480  nd = nj;
481 
482  return ni*Dim;
483  }
484  else
485  {
486  cerr << "NURBSPatch::SetLoopDirection :\n"
487  " Direction error in 2D patch, dir = " << dir << '\n';
488  mfem_error();
489  }
490  }
491  else
492  {
493  if (dir == 0)
494  {
495  sd = Dim;
496  nd = ni;
497 
498  return nj*nk*Dim;
499  }
500  else if (dir == 1)
501  {
502  sd = ni*Dim;
503  nd = nj;
504 
505  return ni*nk*Dim;
506  }
507  else if (dir == 2)
508  {
509  sd = ni*nj*Dim;
510  nd = nk;
511 
512  return ni*nj*Dim;
513  }
514  else
515  {
516  cerr << "NURBSPatch::SetLoopDirection :\n"
517  " Direction error in 3D patch, dir = " << dir << '\n';
518  mfem_error();
519  }
520  }
521 
522  return -1;
523 }
524 
526 {
527  Vector newknots;
528  for (int dir = 0; dir < kv.Size(); dir++)
529  {
530  kv[dir]->UniformRefinement(newknots);
531  KnotInsert(dir, newknots);
532  }
533 }
534 
536 {
537  for (int dir = 0; dir < kv.Size(); dir++)
538  {
539  KnotInsert(dir, *newkv[dir]);
540  }
541 }
542 
543 void NURBSPatch::KnotInsert(int dir, const KnotVector &newkv)
544 {
545  if (dir >= kv.Size() || dir < 0)
546  {
547  mfem_error("NURBSPatch::KnotInsert : Incorrect direction!");
548  }
549 
550  int t = newkv.GetOrder() - kv[dir]->GetOrder();
551 
552  if (t > 0)
553  {
554  DegreeElevate(dir, t);
555  }
556  else if (t < 0)
557  {
558  mfem_error("NURBSPatch::KnotInsert : Incorrect order!");
559  }
560 
561  Vector diff;
562  GetKV(dir)->Difference(newkv, diff);
563  if (diff.Size() > 0)
564  {
565  KnotInsert(dir, diff);
566  }
567 }
568 
569 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
570 void NURBSPatch::KnotInsert(int dir, const Vector &knot)
571 {
572  if (dir >= kv.Size() || dir < 0)
573  {
574  mfem_error("NURBSPatch::KnotInsert : Incorrect direction!");
575  }
576 
577  NURBSPatch &oldp = *this;
578  KnotVector &oldkv = *kv[dir];
579 
580  NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder(),
581  oldkv.GetNCP() + knot.Size());
582  NURBSPatch &newp = *newpatch;
583  KnotVector &newkv = *newp.GetKV(dir);
584 
585  int size = oldp.SetLoopDirection(dir);
586  if (size != newp.SetLoopDirection(dir))
587  {
588  mfem_error("NURBSPatch::KnotInsert : Size mismatch!");
589  }
590 
591  int rr = knot.Size() - 1;
592  int a = oldkv.findKnotSpan(knot(0)) - 1;
593  int b = oldkv.findKnotSpan(knot(rr)) - 1;
594  int pl = oldkv.GetOrder();
595  int ml = oldkv.GetNCP();
596 
597  for (int j = 0; j <= a; j++)
598  {
599  newkv[j] = oldkv[j];
600  }
601  for (int j = b+pl; j <= ml+pl; j++)
602  {
603  newkv[j+rr+1] = oldkv[j];
604  }
605  for (int k = 0; k <= (a-pl); k++)
606  {
607  for (int ll = 0; ll < size; ll++)
608  {
609  newp(k,ll) = oldp(k,ll);
610  }
611  }
612  for (int k = (b-1); k < ml; k++)
613  {
614  for (int ll = 0; ll < size; ll++)
615  {
616  newp(k+rr+1,ll) = oldp(k,ll);
617  }
618  }
619 
620  int i = b+pl-1;
621  int k = b+pl+rr;
622 
623  for (int j = rr; j >= 0; j--)
624  {
625  while ( (knot(j) <= oldkv[i]) && (i > a) )
626  {
627  newkv[k] = oldkv[i];
628  for (int ll = 0; ll < size; ll++)
629  {
630  newp(k-pl-1,ll) = oldp(i-pl-1,ll);
631  }
632 
633  k--;
634  i--;
635  }
636 
637  for (int ll = 0; ll < size; ll++)
638  {
639  newp(k-pl-1,ll) = newp(k-pl,ll);
640  }
641 
642  for (int l = 1; l <= pl; l++)
643  {
644  int ind = k-pl+l;
645  double alfa = newkv[k+l] - knot(j);
646  if (fabs(alfa) == 0.0)
647  {
648  for (int ll = 0; ll < size; ll++)
649  {
650  newp(ind-1,ll) = newp(ind,ll);
651  }
652  }
653  else
654  {
655  alfa = alfa/(newkv[k+l] - oldkv[i-pl+l]);
656  for (int ll = 0; ll < size; ll++)
657  {
658  newp(ind-1,ll) = alfa*newp(ind-1,ll) + (1.0-alfa)*newp(ind,ll);
659  }
660  }
661  }
662 
663  newkv[k] = knot(j);
664  k--;
665  }
666 
667  newkv.GetElements();
668 
669  swap(newpatch);
670 }
671 
673 {
674  for (int dir = 0; dir < kv.Size(); dir++)
675  {
676  DegreeElevate(dir, t);
677  }
678 }
679 
680 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
681 void NURBSPatch::DegreeElevate(int dir, int t)
682 {
683  if (dir >= kv.Size() || dir < 0)
684  {
685  mfem_error("NURBSPatch::DegreeElevate : Incorrect direction!");
686  }
687 
688  int i, j, k, kj, mpi, mul, mh, kind, cind, first, last;
689  int r, a, b, oldr, save, s, tr, lbz, rbz, l;
690  double inv, ua, ub, numer, alf, den, bet, gam;
691 
692  NURBSPatch &oldp = *this;
693  KnotVector &oldkv = *kv[dir];
694 
695  NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder() + t,
696  oldkv.GetNCP() + oldkv.GetNE()*t);
697  NURBSPatch &newp = *newpatch;
698  KnotVector &newkv = *newp.GetKV(dir);
699 
700  int size = oldp.SetLoopDirection(dir);
701  if (size != newp.SetLoopDirection(dir))
702  {
703  mfem_error("NURBSPatch::DegreeElevate : Size mismatch!");
704  }
705 
706  int p = oldkv.GetOrder();
707  int n = oldkv.GetNCP()-1;
708 
709  DenseMatrix bezalfs (p+t+1, p+1);
710  DenseMatrix bpts (p+1, size);
711  DenseMatrix ebpts (p+t+1, size);
712  DenseMatrix nextbpts(p-1, size);
713  Vector alphas (p-1);
714 
715  int m = n + p + 1;
716  int ph = p + t;
717  int ph2 = ph/2;
718 
719  {
720  Array2D<int> binom(ph+1, ph+1);
721  for (i = 0; i <= ph; i++)
722  {
723  binom(i,0) = binom(i,i) = 1;
724  for (j = 1; j < i; j++)
725  {
726  binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
727  }
728  }
729 
730  bezalfs(0,0) = 1.0;
731  bezalfs(ph,p) = 1.0;
732 
733  for (i = 1; i <= ph2; i++)
734  {
735  inv = 1.0/binom(ph,i);
736  mpi = min(p,i);
737  for (j = max(0,i-t); j <= mpi; j++)
738  {
739  bezalfs(i,j) = inv*binom(p,j)*binom(t,i-j);
740  }
741  }
742  }
743 
744  for (i = ph2+1; i < ph; i++)
745  {
746  mpi = min(p,i);
747  for (j = max(0,i-t); j <= mpi; j++)
748  {
749  bezalfs(i,j) = bezalfs(ph-i,p-j);
750  }
751  }
752 
753  mh = ph;
754  kind = ph + 1;
755  r = -1;
756  a = p;
757  b = p + 1;
758  cind = 1;
759  ua = oldkv[0];
760  for (l = 0; l < size; l++)
761  {
762  newp(0,l) = oldp(0,l);
763  }
764  for (i = 0; i <= ph; i++)
765  {
766  newkv[i] = ua;
767  }
768 
769  for (i = 0; i <= p; i++)
770  {
771  for (l = 0; l < size; l++)
772  {
773  bpts(i,l) = oldp(i,l);
774  }
775  }
776 
777  while (b < m)
778  {
779  i = b;
780  while (b < m && oldkv[b] == oldkv[b+1]) { b++; }
781 
782  mul = b-i+1;
783 
784  mh = mh + mul + t;
785  ub = oldkv[b];
786  oldr = r;
787  r = p-mul;
788  if (oldr > 0) { lbz = (oldr+2)/2; }
789  else { lbz = 1; }
790 
791  if (r > 0) { rbz = ph-(r+1)/2; }
792  else { rbz = ph; }
793 
794  if (r > 0)
795  {
796  numer = ub - ua;
797  for (k = p ; k > mul; k--)
798  {
799  alphas[k-mul-1] = numer/(oldkv[a+k]-ua);
800  }
801 
802  for (j = 1; j <= r; j++)
803  {
804  save = r-j;
805  s = mul+j;
806  for (k = p; k >= s; k--)
807  {
808  for (l = 0; l < size; l++)
809  bpts(k,l) = (alphas[k-s]*bpts(k,l) +
810  (1.0-alphas[k-s])*bpts(k-1,l));
811  }
812  for (l = 0; l < size; l++)
813  {
814  nextbpts(save,l) = bpts(p,l);
815  }
816  }
817  }
818 
819  for (i = lbz; i <= ph; i++)
820  {
821  for (l = 0; l < size; l++)
822  {
823  ebpts(i,l) = 0.0;
824  }
825  mpi = min(p,i);
826  for (j = max(0,i-t); j <= mpi; j++)
827  {
828  for (l = 0; l < size; l++)
829  {
830  ebpts(i,l) += bezalfs(i,j)*bpts(j,l);
831  }
832  }
833  }
834 
835  if (oldr > 1)
836  {
837  first = kind-2;
838  last = kind;
839  den = ub-ua;
840  bet = (ub-newkv[kind-1])/den;
841 
842  for (tr = 1; tr < oldr; tr++)
843  {
844  i = first;
845  j = last;
846  kj = j-kind+1;
847  while (j-i > tr)
848  {
849  if (i < cind)
850  {
851  alf = (ub-newkv[i])/(ua-newkv[i]);
852  for (l = 0; l < size; l++)
853  {
854  newp(i,l) = alf*newp(i,l)-(1.0-alf)*newp(i-1,l);
855  }
856  }
857  if (j >= lbz)
858  {
859  if ((j-tr) <= (kind-ph+oldr))
860  {
861  gam = (ub-newkv[j-tr])/den;
862  for (l = 0; l < size; l++)
863  {
864  ebpts(kj,l) = gam*ebpts(kj,l) + (1.0-gam)*ebpts(kj+1,l);
865  }
866  }
867  else
868  {
869  for (l = 0; l < size; l++)
870  {
871  ebpts(kj,l) = bet*ebpts(kj,l) + (1.0-bet)*ebpts(kj+1,l);
872  }
873  }
874  }
875  i = i+1;
876  j = j-1;
877  kj = kj-1;
878  }
879  first--;
880  last++;
881  }
882  }
883 
884  if (a != p)
885  {
886  for (i = 0; i < (ph-oldr); i++)
887  {
888  newkv[kind] = ua;
889  kind = kind+1;
890  }
891  }
892  for (j = lbz; j <= rbz; j++)
893  {
894  for (l = 0; l < size; l++)
895  {
896  newp(cind,l) = ebpts(j,l);
897  }
898  cind = cind +1;
899  }
900 
901  if (b < m)
902  {
903  for (j = 0; j <r; j++)
904  for (l = 0; l < size; l++)
905  {
906  bpts(j,l) = nextbpts(j,l);
907  }
908 
909  for (j = r; j <= p; j++)
910  for (l = 0; l < size; l++)
911  {
912  bpts(j,l) = oldp(b-p+j,l);
913  }
914 
915  a = b;
916  b = b+1;
917  ua = ub;
918  }
919  else
920  {
921  for (i = 0; i <= ph; i++)
922  {
923  newkv[kind+i] = ub;
924  }
925  }
926  }
927  newkv.GetElements();
928 
929  swap(newpatch);
930 }
931 
933 {
934  int size = SetLoopDirection(dir);
935 
936  for (int id = 0; id < nd/2; id++)
937  for (int i = 0; i < size; i++)
938  {
939  Swap<double>((*this)(id,i), (*this)(nd-1-id,i));
940  }
941  kv[dir]->Flip();
942 }
943 
944 void NURBSPatch::SwapDirections(int dir1, int dir2)
945 {
946  if (abs(dir1-dir2) == 2)
947  mfem_error("NURBSPatch::SwapDirections :"
948  " directions 0 and 2 are not supported!");
949 
951 
952  kv.Copy(nkv);
953  Swap<KnotVector *>(nkv[dir1], nkv[dir2]);
954  NURBSPatch *newpatch = new NURBSPatch(nkv, Dim);
955 
956  int size = SetLoopDirection(dir1);
957  newpatch->SetLoopDirection(dir2);
958 
959  for (int id = 0; id < nd; id++)
960  for (int i = 0; i < size; i++)
961  {
962  (*newpatch)(id,i) = (*this)(id,i);
963  }
964 
965  swap(newpatch);
966 }
967 
968 void NURBSPatch::Get3DRotationMatrix(double n[], double angle, double r,
969  DenseMatrix &T)
970 {
971  double c, s, c1;
972  double l2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
973  double l = sqrt(l2);
974 
975  if (fabs(angle) == M_PI_2)
976  {
977  s = r*copysign(1., angle);
978  c = 0.;
979  c1 = -1.;
980  }
981  else if (fabs(angle) == M_PI)
982  {
983  s = 0.;
984  c = -r;
985  c1 = c - 1.;
986  }
987  else
988  {
989  s = r*sin(angle);
990  c = r*cos(angle);
991  c1 = c - 1.;
992  }
993 
994  T.SetSize(3);
995 
996  T(0,0) = (n[0]*n[0] + (n[1]*n[1] + n[2]*n[2])*c)/l2;
997  T(0,1) = -(n[0]*n[1]*c1)/l2 - (n[2]*s)/l;
998  T(0,2) = -(n[0]*n[2]*c1)/l2 + (n[1]*s)/l;
999  T(1,0) = -(n[0]*n[1]*c1)/l2 + (n[2]*s)/l;
1000  T(1,1) = (n[1]*n[1] + (n[0]*n[0] + n[2]*n[2])*c)/l2;
1001  T(1,2) = -(n[1]*n[2]*c1)/l2 - (n[0]*s)/l;
1002  T(2,0) = -(n[0]*n[2]*c1)/l2 - (n[1]*s)/l;
1003  T(2,1) = -(n[1]*n[2]*c1)/l2 + (n[0]*s)/l;
1004  T(2,2) = (n[2]*n[2] + (n[0]*n[0] + n[1]*n[1])*c)/l2;
1005 }
1006 
1007 void NURBSPatch::Rotate3D(double n[], double angle)
1008 {
1009  if (Dim != 4)
1010  {
1011  mfem_error("NURBSPatch::Rotate3D : not a NURBSPatch in 3D!");
1012  }
1013 
1014  DenseMatrix T(3);
1015  Vector x(3), y(NULL, 3);
1016 
1017  Get3DRotationMatrix(n, angle, 1., T);
1018 
1019  int size = 1;
1020  for (int i = 0; i < kv.Size(); i++)
1021  {
1022  size *= kv[i]->GetNCP();
1023  }
1024 
1025  for (int i = 0; i < size; i++)
1026  {
1027  y.SetData(data + i*Dim);
1028  x = y;
1029  T.Mult(x, y);
1030  }
1031 }
1032 
1034 {
1035  int maxd = -1;
1036 
1037  for (int dir = 0; dir < kv.Size(); dir++)
1038  if (maxd < kv[dir]->GetOrder())
1039  {
1040  maxd = kv[dir]->GetOrder();
1041  }
1042 
1043  for (int dir = 0; dir < kv.Size(); dir++)
1044  if (maxd > kv[dir]->GetOrder())
1045  {
1046  DegreeElevate(dir, maxd - kv[dir]->GetOrder());
1047  }
1048 
1049  return maxd;
1050 }
1051 
1053 {
1054  if (p1.kv.Size() != p2.kv.Size() || p1.Dim != p2.Dim)
1055  {
1056  mfem_error("Interpolate(NURBSPatch &, NURBSPatch &)");
1057  }
1058 
1059  int size = 1, dim = p1.Dim;
1060  Array<KnotVector *> kv(p1.kv.Size() + 1);
1061 
1062  for (int i = 0; i < p1.kv.Size(); i++)
1063  {
1064  if (p1.kv[i]->GetOrder() < p2.kv[i]->GetOrder())
1065  {
1066  p1.KnotInsert(i, *p2.kv[i]);
1067  p2.KnotInsert(i, *p1.kv[i]);
1068  }
1069  else
1070  {
1071  p2.KnotInsert(i, *p1.kv[i]);
1072  p1.KnotInsert(i, *p2.kv[i]);
1073  }
1074  kv[i] = p1.kv[i];
1075  size *= kv[i]->GetNCP();
1076  }
1077 
1078  kv.Last() = new KnotVector(1, 2);
1079  KnotVector &nkv = *kv.Last();
1080  nkv[0] = nkv[1] = 0.0;
1081  nkv[2] = nkv[3] = 1.0;
1082  nkv.GetElements();
1083 
1084  NURBSPatch *patch = new NURBSPatch(kv, dim);
1085  delete kv.Last();
1086 
1087  for (int i = 0; i < size; i++)
1088  {
1089  for (int d = 0; d < dim; d++)
1090  {
1091  patch->data[i*dim+d] = p1.data[i*dim+d];
1092  patch->data[(i+size)*dim+d] = p2.data[i*dim+d];
1093  }
1094  }
1095 
1096  return patch;
1097 }
1098 
1099 NURBSPatch *Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
1100 {
1101  if (patch.Dim != 4)
1102  {
1103  mfem_error("Revolve3D(NURBSPatch &, double [], double)");
1104  }
1105 
1106  int size = 1, ns;
1107  Array<KnotVector *> nkv(patch.kv.Size() + 1);
1108 
1109  for (int i = 0; i < patch.kv.Size(); i++)
1110  {
1111  nkv[i] = patch.kv[i];
1112  size *= nkv[i]->GetNCP();
1113  }
1114  ns = 2*times + 1;
1115  nkv.Last() = new KnotVector(2, ns);
1116  KnotVector &lkv = *nkv.Last();
1117  lkv[0] = lkv[1] = lkv[2] = 0.0;
1118  for (int i = 1; i < times; i++)
1119  {
1120  lkv[2*i+1] = lkv[2*i+2] = i;
1121  }
1122  lkv[ns] = lkv[ns+1] = lkv[ns+2] = times;
1123  lkv.GetElements();
1124  NURBSPatch *newpatch = new NURBSPatch(nkv, 4);
1125  delete nkv.Last();
1126 
1127  DenseMatrix T(3), T2(3);
1128  Vector u(NULL, 3), v(NULL, 3);
1129 
1130  NURBSPatch::Get3DRotationMatrix(n, ang, 1., T);
1131  double c = cos(ang/2);
1132  NURBSPatch::Get3DRotationMatrix(n, ang/2, 1./c, T2);
1133  T2 *= c;
1134 
1135  double *op = patch.data, *np;
1136  for (int i = 0; i < size; i++)
1137  {
1138  np = newpatch->data + 4*i;
1139  for (int j = 0; j < 4; j++)
1140  {
1141  np[j] = op[j];
1142  }
1143  for (int j = 0; j < times; j++)
1144  {
1145  u.SetData(np);
1146  v.SetData(np += 4*size);
1147  T2.Mult(u, v);
1148  v[3] = c*u[3];
1149  v.SetData(np += 4*size);
1150  T.Mult(u, v);
1151  v[3] = u[3];
1152  }
1153  op += 4;
1154  }
1155 
1156  return newpatch;
1157 }
1158 
1159 
1161 {
1162  // Read topology
1163  patchTopo = new Mesh;
1164  patchTopo->LoadPatchTopo(input, edge_to_knot);
1165  own_topo = 1;
1166 
1167  CheckPatches();
1168  // CheckBdrPatches();
1169 
1170  Mesh::skip_comment_lines(input, '#');
1171 
1172  // Read knotvectors or patches
1173  string ident;
1174  input >> ws >> ident; // 'knotvectors' or 'patches'
1175  if (ident == "knotvectors")
1176  {
1177  input >> NumOfKnotVectors;
1178  knotVectors.SetSize(NumOfKnotVectors);
1179  for (int i = 0; i < NumOfKnotVectors; i++)
1180  {
1181  knotVectors[i] = new KnotVector(input);
1182  if (knotVectors[i]->GetOrder() != knotVectors[0]->GetOrder())
1183  mfem_error("NURBSExtension::NURBSExtension :\n"
1184  " Variable orders are not supported!");
1185  }
1186  Order = knotVectors[0]->GetOrder();
1187  }
1188  else if (ident == "patches")
1189  {
1190  patches.SetSize(GetNP());
1191  for (int p = 0; p < patches.Size(); p++)
1192  {
1193  Mesh::skip_comment_lines(input, '#');
1194  patches[p] = new NURBSPatch(input);
1195  }
1196 
1197  NumOfKnotVectors = 0;
1198  for (int i = 0; i < patchTopo->GetNEdges(); i++)
1199  if (NumOfKnotVectors < KnotInd(i))
1200  {
1201  NumOfKnotVectors = KnotInd(i);
1202  }
1203  NumOfKnotVectors++;
1204  knotVectors.SetSize(NumOfKnotVectors);
1205  knotVectors = NULL;
1206 
1207  Array<int> edges, oedge;
1208  for (int p = 0; p < patches.Size(); p++)
1209  {
1210  patchTopo->GetElementEdges(p, edges, oedge);
1211  if (Dimension() == 2)
1212  {
1213  if (knotVectors[KnotInd(edges[0])] == NULL)
1214  knotVectors[KnotInd(edges[0])] =
1215  new KnotVector(*patches[p]->GetKV(0));
1216  if (knotVectors[KnotInd(edges[1])] == NULL)
1217  knotVectors[KnotInd(edges[1])] =
1218  new KnotVector(*patches[p]->GetKV(1));
1219  }
1220  else
1221  {
1222  if (knotVectors[KnotInd(edges[0])] == NULL)
1223  knotVectors[KnotInd(edges[0])] =
1224  new KnotVector(*patches[p]->GetKV(0));
1225  if (knotVectors[KnotInd(edges[3])] == NULL)
1226  knotVectors[KnotInd(edges[3])] =
1227  new KnotVector(*patches[p]->GetKV(1));
1228  if (knotVectors[KnotInd(edges[8])] == NULL)
1229  knotVectors[KnotInd(edges[8])] =
1230  new KnotVector(*patches[p]->GetKV(2));
1231  }
1232  }
1233  Order = knotVectors[0]->GetOrder();
1234  }
1235 
1236  GenerateOffsets();
1237  CountElements();
1238  CountBdrElements();
1239  // NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs
1240 
1241  Mesh::skip_comment_lines(input, '#');
1242 
1243  // Check for a list of mesh elements
1244  if (patches.Size() == 0)
1245  {
1246  input >> ws >> ident;
1247  }
1248  if (patches.Size() == 0 && ident == "mesh_elements")
1249  {
1250  input >> NumOfActiveElems;
1251  activeElem.SetSize(GetGNE());
1252  activeElem = false;
1253  int glob_elem;
1254  for (int i = 0; i < NumOfActiveElems; i++)
1255  {
1256  input >> glob_elem;
1257  activeElem[glob_elem] = true;
1258  }
1259 
1260  Mesh::skip_comment_lines(input, '#');
1261  input >> ws >> ident;
1262  }
1263  else
1264  {
1265  NumOfActiveElems = NumOfElements;
1266  activeElem.SetSize(NumOfElements);
1267  activeElem = true;
1268  }
1269 
1270  GenerateActiveVertices();
1271  GenerateElementDofTable();
1272  GenerateActiveBdrElems();
1273  GenerateBdrElementDofTable();
1274 
1275  if (patches.Size() == 0)
1276  {
1277  // weights
1278  if (ident == "weights")
1279  {
1280  weights.Load(input, GetNDof());
1281  }
1282  else // e.g. ident = "unitweights" or "autoweights"
1283  {
1284  weights.SetSize(GetNDof());
1285  weights = 1.0;
1286  }
1287  }
1288 }
1289 
1291 {
1292  Order = Order_;
1293 
1294  patchTopo = parent->patchTopo;
1295  own_topo = 0;
1296 
1297  parent->edge_to_knot.Copy(edge_to_knot);
1298 
1299  NumOfKnotVectors = parent->GetNKV();
1300  knotVectors.SetSize(NumOfKnotVectors);
1301  for (int i = 0; i < NumOfKnotVectors; i++)
1302  {
1303  knotVectors[i] =
1304  parent->GetKnotVector(i)->DegreeElevate(Order-parent->GetOrder());
1305  }
1306 
1307  // copy some data from parent
1308  NumOfElements = parent->NumOfElements;
1309  NumOfBdrElements = parent->NumOfBdrElements;
1310 
1311  GenerateOffsets(); // dof offsets will be different from parent
1312 
1313  NumOfActiveVertices = parent->NumOfActiveVertices;
1314  NumOfActiveElems = parent->NumOfActiveElems;
1315  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1316  parent->activeVert.Copy(activeVert);
1317  parent->activeElem.Copy(activeElem);
1318  parent->activeBdrElem.Copy(activeBdrElem);
1319 
1320  GenerateElementDofTable();
1321  GenerateBdrElementDofTable();
1322 
1323  weights.SetSize(GetNDof());
1324  weights = 1.0;
1325 }
1326 
1327 NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
1328 {
1329  NURBSExtension *parent = mesh_array[0]->NURBSext;
1330 
1331  if (!parent->own_topo)
1332  mfem_error("NURBSExtension::NURBSExtension :\n"
1333  " parent does not own the patch topology!");
1334  patchTopo = parent->patchTopo;
1335  own_topo = 1;
1336  parent->own_topo = 0;
1337 
1338  parent->edge_to_knot.Copy(edge_to_knot);
1339 
1340  Order = parent->GetOrder();
1341 
1342  NumOfKnotVectors = parent->GetNKV();
1343  knotVectors.SetSize(NumOfKnotVectors);
1344  for (int i = 0; i < NumOfKnotVectors; i++)
1345  {
1346  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1347  }
1348 
1349  GenerateOffsets();
1350  CountElements();
1351  CountBdrElements();
1352 
1353  // assuming the meshes define a partitioning of all the elements
1354  NumOfActiveElems = NumOfElements;
1355  activeElem.SetSize(NumOfElements);
1356  activeElem = true;
1357 
1358  GenerateActiveVertices();
1359  GenerateElementDofTable();
1360  GenerateActiveBdrElems();
1361  GenerateBdrElementDofTable();
1362 
1363  weights.SetSize(GetNDof());
1364  MergeWeights(mesh_array, num_pieces);
1365 }
1366 
1368 {
1369  if (patches.Size() == 0)
1370  {
1371  delete bel_dof;
1372  delete el_dof;
1373  }
1374 
1375  for (int i = 0; i < knotVectors.Size(); i++)
1376  {
1377  delete knotVectors[i];
1378  }
1379 
1380  for (int i = 0; i < patches.Size(); i++)
1381  {
1382  delete patches[i];
1383  }
1384 
1385  if (own_topo)
1386  {
1387  delete patchTopo;
1388  }
1389 }
1390 
1391 void NURBSExtension::Print(std::ostream &out) const
1392 {
1393  patchTopo->PrintTopo(out, edge_to_knot);
1394  if (patches.Size() == 0)
1395  {
1396  out << "\nknotvectors\n" << NumOfKnotVectors << '\n';
1397  for (int i = 0; i < NumOfKnotVectors; i++)
1398  {
1399  knotVectors[i]->Print(out);
1400  }
1401 
1402  if (NumOfActiveElems < NumOfElements)
1403  {
1404  out << "\nmesh_elements\n" << NumOfActiveElems << '\n';
1405  for (int i = 0; i < NumOfElements; i++)
1406  if (activeElem[i])
1407  {
1408  out << i << '\n';
1409  }
1410  }
1411 
1412  out << "\nweights\n";
1413  weights.Print(out, 1);
1414  }
1415  else
1416  {
1417  out << "\npatches\n";
1418  for (int p = 0; p < patches.Size(); p++)
1419  {
1420  out << "\n# patch " << p << "\n\n";
1421  patches[p]->Print(out);
1422  }
1423  }
1424 }
1425 
1427 {
1428  out <<
1429  "NURBS Mesh entity sizes:\n"
1430  "Dimension = " << Dimension() << "\n"
1431  "Order = " << GetOrder() << "\n"
1432  "NumOfKnotVectors = " << GetNKV() << "\n"
1433  "NumOfPatches = " << GetNP() << "\n"
1434  "NumOfBdrPatches = " << GetNBP() << "\n"
1435  "NumOfVertices = " << GetGNV() << "\n"
1436  "NumOfElements = " << GetGNE() << "\n"
1437  "NumOfBdrElements = " << GetGNBE() << "\n"
1438  "NumOfDofs = " << GetNTotalDof() << "\n"
1439  "NumOfActiveVertices = " << GetNV() << "\n"
1440  "NumOfActiveElems = " << GetNE() << "\n"
1441  "NumOfActiveBdrElems = " << GetNBE() << "\n"
1442  "NumOfActiveDofs = " << GetNDof() << '\n';
1443  for (int i = 0; i < NumOfKnotVectors; i++)
1444  {
1445  out << ' ' << i + 1 << ") ";
1446  knotVectors[i]->Print(out);
1447  }
1448  cout << endl;
1449 }
1450 
1452 {
1453  int vert[8], nv, g_el, nx, ny, nz, dim = Dimension();
1454 
1455  NURBSPatchMap p2g(this);
1456  KnotVector *kv[3];
1457 
1458  g_el = 0;
1459  activeVert.SetSize(GetGNV());
1460  activeVert = -1;
1461  for (int p = 0; p < GetNP(); p++)
1462  {
1463  p2g.SetPatchVertexMap(p, kv);
1464 
1465  nx = p2g.nx();
1466  ny = p2g.ny();
1467  nz = (dim == 3) ? p2g.nz() : 1;
1468 
1469  for (int k = 0; k < nz; k++)
1470  {
1471  for (int j = 0; j < ny; j++)
1472  {
1473  for (int i = 0; i < nx; i++)
1474  {
1475  if (activeElem[g_el])
1476  {
1477  if (dim == 2)
1478  {
1479  vert[0] = p2g(i, j );
1480  vert[1] = p2g(i+1,j );
1481  vert[2] = p2g(i+1,j+1);
1482  vert[3] = p2g(i, j+1);
1483  nv = 4;
1484  }
1485  else
1486  {
1487  vert[0] = p2g(i, j, k);
1488  vert[1] = p2g(i+1,j, k);
1489  vert[2] = p2g(i+1,j+1,k);
1490  vert[3] = p2g(i, j+1,k);
1491 
1492  vert[4] = p2g(i, j, k+1);
1493  vert[5] = p2g(i+1,j, k+1);
1494  vert[6] = p2g(i+1,j+1,k+1);
1495  vert[7] = p2g(i, j+1,k+1);
1496  nv = 8;
1497  }
1498 
1499  for (int v = 0; v < nv; v++)
1500  {
1501  activeVert[vert[v]] = 1;
1502  }
1503  }
1504  g_el++;
1505  }
1506  }
1507  }
1508  }
1509 
1510  NumOfActiveVertices = 0;
1511  for (int i = 0; i < GetGNV(); i++)
1512  if (activeVert[i] == 1)
1513  {
1514  activeVert[i] = NumOfActiveVertices++;
1515  }
1516 }
1517 
1519 {
1520  int dim = Dimension();
1521  Array<KnotVector *> kv(dim);
1522 
1523  activeBdrElem.SetSize(GetGNBE());
1524  if (GetGNE() == GetNE())
1525  {
1526  activeBdrElem = true;
1527  NumOfActiveBdrElems = GetGNBE();
1528  return;
1529  }
1530  activeBdrElem = false;
1531  NumOfActiveBdrElems = 0;
1532  // the mesh will generate the actual boundary including boundary
1533  // elements that are not on boundary patches. we use this for
1534  // visualization of processor boundaries
1535 
1536  // TODO: generate actual boundary?
1537 }
1538 
1539 
1540 void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces)
1541 {
1542  Array<int> lelem_elem;
1543 
1544  for (int i = 0; i < num_pieces; i++)
1545  {
1546  NURBSExtension *lext = mesh_array[i]->NURBSext;
1547 
1548  lext->GetElementLocalToGlobal(lelem_elem);
1549 
1550  for (int lel = 0; lel < lext->GetNE(); lel++)
1551  {
1552  int gel = lelem_elem[lel];
1553 
1554  int nd = el_dof->RowSize(gel);
1555  int *gdofs = el_dof->GetRow(gel);
1556  int *ldofs = lext->el_dof->GetRow(lel);
1557  for (int j = 0; j < nd; j++)
1558  {
1559  weights(gdofs[j]) = lext->weights(ldofs[j]);
1560  }
1561  }
1562  }
1563 }
1564 
1566  GridFunction *gf_array[], int num_pieces, GridFunction &merged)
1567 {
1568  FiniteElementSpace *gfes = merged.FESpace();
1569  Array<int> lelem_elem, dofs;
1570  Vector lvec;
1571 
1572  for (int i = 0; i < num_pieces; i++)
1573  {
1574  FiniteElementSpace *lfes = gf_array[i]->FESpace();
1575  NURBSExtension *lext = lfes->GetMesh()->NURBSext;
1576 
1577  lext->GetElementLocalToGlobal(lelem_elem);
1578 
1579  for (int lel = 0; lel < lext->GetNE(); lel++)
1580  {
1581  lfes->GetElementVDofs(lel, dofs);
1582  gf_array[i]->GetSubVector(dofs, lvec);
1583 
1584  gfes->GetElementVDofs(lelem_elem[lel], dofs);
1585  merged.SetSubVector(dofs, lvec);
1586  }
1587  }
1588 }
1589 
1591 {
1592  Array<int> edges;
1593  Array<int> oedge;
1594 
1595  for (int p = 0; p < GetNP(); p++)
1596  {
1597  patchTopo->GetElementEdges(p, edges, oedge);
1598 
1599  for (int i = 0; i < edges.Size(); i++)
1600  {
1601  edges[i] = edge_to_knot[edges[i]];
1602  if (oedge[i] < 0)
1603  {
1604  edges[i] = -1 - edges[i];
1605  }
1606  }
1607 
1608  if ((Dimension() == 2 &&
1609  (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
1610 
1611  (Dimension() == 3 &&
1612  (edges[0] != edges[2] || edges[0] != edges[4] ||
1613  edges[0] != edges[6] || edges[1] != edges[3] ||
1614  edges[1] != edges[5] || edges[1] != edges[7] ||
1615  edges[8] != edges[9] || edges[8] != edges[10] ||
1616  edges[8] != edges[11])))
1617  {
1618  cerr << "NURBSExtension::CheckPatch (patch = " << p
1619  << ")\n Inconsistent edge-to-knot mapping!\n";
1620  mfem_error();
1621  }
1622 
1623  if ((Dimension() == 2 &&
1624  (edges[0] < 0 || edges[1] < 0)) ||
1625 
1626  (Dimension() == 3 &&
1627  (edges[0] < 0 || edges[3] < 0 || edges[8] < 0)))
1628  {
1629  cerr << "NURBSExtension::CheckPatch (patch = " << p
1630  << ") : Bad orientation!\n";
1631  mfem_error();
1632  }
1633  }
1634 }
1635 
1637 {
1638  Array<int> edges;
1639  Array<int> oedge;
1640 
1641  for (int p = 0; p < GetNBP(); p++)
1642  {
1643  patchTopo->GetBdrElementEdges(p, edges, oedge);
1644 
1645  for (int i = 0; i < edges.Size(); i++)
1646  {
1647  edges[i] = edge_to_knot[edges[i]];
1648  if (oedge[i] < 0)
1649  {
1650  edges[i] = -1 - edges[i];
1651  }
1652  }
1653 
1654  if ((Dimension() == 2 && (edges[0] < 0)) ||
1655  (Dimension() == 3 && (edges[0] < 0 || edges[1] < 0)))
1656  {
1657  cerr << "NURBSExtension::CheckBdrPatch (boundary patch = "
1658  << p << ") : Bad orientation!\n";
1659  mfem_error();
1660  }
1661  }
1662 }
1663 
1665 {
1666  Array<int> edges;
1667  Array<int> orient;
1668 
1669  kv.SetSize(Dimension());
1670  patchTopo->GetElementEdges(p, edges, orient);
1671 
1672  if (Dimension() == 2)
1673  {
1674  kv[0] = KnotVec(edges[0]);
1675  kv[1] = KnotVec(edges[1]);
1676  }
1677  else
1678  {
1679  kv[0] = KnotVec(edges[0]);
1680  kv[1] = KnotVec(edges[3]);
1681  kv[2] = KnotVec(edges[8]);
1682  }
1683 }
1684 
1686 {
1687  Array<int> edges;
1688  Array<int> orient;
1689 
1690  kv.SetSize(Dimension() - 1);
1691  patchTopo->GetBdrElementEdges(p, edges, orient);
1692 
1693  if (Dimension() == 2)
1694  {
1695  kv[0] = KnotVec(edges[0]);
1696  }
1697  else
1698  {
1699  kv[0] = KnotVec(edges[0]);
1700  kv[1] = KnotVec(edges[1]);
1701  }
1702 }
1703 
1705 {
1706  int nv = patchTopo->GetNV();
1707  int ne = patchTopo->GetNEdges();
1708  int nf = patchTopo->GetNFaces();
1709  int np = patchTopo->GetNE();
1710  int meshCounter, spaceCounter, dim = Dimension();
1711 
1712  Array<int> edges;
1713  Array<int> orient;
1714 
1715  v_meshOffsets.SetSize(nv);
1716  e_meshOffsets.SetSize(ne);
1717  f_meshOffsets.SetSize(nf);
1718  p_meshOffsets.SetSize(np);
1719 
1720  v_spaceOffsets.SetSize(nv);
1721  e_spaceOffsets.SetSize(ne);
1722  f_spaceOffsets.SetSize(nf);
1723  p_spaceOffsets.SetSize(np);
1724 
1725  // Get vertex offsets
1726  for (meshCounter = 0; meshCounter < nv; meshCounter++)
1727  {
1728  v_meshOffsets[meshCounter] = meshCounter;
1729  v_spaceOffsets[meshCounter] = meshCounter;
1730  }
1731  spaceCounter = meshCounter;
1732 
1733  // Get edge offsets
1734  for (int e = 0; e < ne; e++)
1735  {
1736  e_meshOffsets[e] = meshCounter;
1737  e_spaceOffsets[e] = spaceCounter;
1738  meshCounter += KnotVec(e)->GetNE() - 1;
1739  spaceCounter += KnotVec(e)->GetNCP() - 2;
1740  }
1741 
1742  // Get face offsets
1743  for (int f = 0; f < nf; f++)
1744  {
1745  f_meshOffsets[f] = meshCounter;
1746  f_spaceOffsets[f] = spaceCounter;
1747 
1748  patchTopo->GetFaceEdges(f, edges, orient);
1749 
1750  meshCounter +=
1751  (KnotVec(edges[0])->GetNE() - 1) *
1752  (KnotVec(edges[1])->GetNE() - 1);
1753  spaceCounter +=
1754  (KnotVec(edges[0])->GetNCP() - 2) *
1755  (KnotVec(edges[1])->GetNCP() - 2);
1756  }
1757 
1758  // Get patch offsets
1759  for (int p = 0; p < np; p++)
1760  {
1761  p_meshOffsets[p] = meshCounter;
1762  p_spaceOffsets[p] = spaceCounter;
1763 
1764  patchTopo->GetElementEdges(p, edges, orient);
1765 
1766  if (dim == 2)
1767  {
1768  meshCounter +=
1769  (KnotVec(edges[0])->GetNE() - 1) *
1770  (KnotVec(edges[1])->GetNE() - 1);
1771  spaceCounter +=
1772  (KnotVec(edges[0])->GetNCP() - 2) *
1773  (KnotVec(edges[1])->GetNCP() - 2);
1774  }
1775  else
1776  {
1777  meshCounter +=
1778  (KnotVec(edges[0])->GetNE() - 1) *
1779  (KnotVec(edges[3])->GetNE() - 1) *
1780  (KnotVec(edges[8])->GetNE() - 1);
1781  spaceCounter +=
1782  (KnotVec(edges[0])->GetNCP() - 2) *
1783  (KnotVec(edges[3])->GetNCP() - 2) *
1784  (KnotVec(edges[8])->GetNCP() - 2);
1785  }
1786  }
1787  NumOfVertices = meshCounter;
1788  NumOfDofs = spaceCounter;
1789 }
1790 
1792 {
1793  int dim = Dimension();
1794  Array<KnotVector *> kv(dim);
1795 
1796  NumOfElements = 0;
1797  for (int p = 0; p < GetNP(); p++)
1798  {
1799  GetPatchKnotVectors(p, kv);
1800 
1801  int ne = kv[0]->GetNE();
1802  for (int d = 1; d < dim; d++)
1803  {
1804  ne *= kv[d]->GetNE();
1805  }
1806 
1807  NumOfElements += ne;
1808  }
1809 }
1810 
1812 {
1813  int dim = Dimension() - 1;
1814  Array<KnotVector *> kv(dim);
1815 
1816  NumOfBdrElements = 0;
1817  for (int p = 0; p < GetNBP(); p++)
1818  {
1819  GetBdrPatchKnotVectors(p, kv);
1820 
1821  int ne = kv[0]->GetNE();
1822  for (int d = 1; d < dim; d++)
1823  {
1824  ne *= kv[d]->GetNE();
1825  }
1826 
1827  NumOfBdrElements += ne;
1828  }
1829 }
1830 
1832 {
1833  elements.SetSize(GetNE());
1834 
1835  if (Dimension() == 2)
1836  {
1837  Get2DElementTopo(elements);
1838  }
1839  else
1840  {
1841  Get3DElementTopo(elements);
1842  }
1843 }
1844 
1846 {
1847  int el = 0;
1848  int eg = 0;
1849  int ind[4];
1850  NURBSPatchMap p2g(this);
1851  KnotVector *kv[2];
1852 
1853  for (int p = 0; p < GetNP(); p++)
1854  {
1855  p2g.SetPatchVertexMap(p, kv);
1856  int nx = p2g.nx();
1857  int ny = p2g.ny();
1858 
1859  int patch_attr = patchTopo->GetAttribute(p);
1860 
1861  for (int j = 0; j < ny; j++)
1862  {
1863  for (int i = 0; i < nx; i++)
1864  {
1865  if (activeElem[eg])
1866  {
1867  ind[0] = activeVert[p2g(i, j )];
1868  ind[1] = activeVert[p2g(i+1,j )];
1869  ind[2] = activeVert[p2g(i+1,j+1)];
1870  ind[3] = activeVert[p2g(i, j+1)];
1871 
1872  elements[el] = new Quadrilateral(ind, patch_attr);
1873  el++;
1874  }
1875  eg++;
1876  }
1877  }
1878  }
1879 }
1880 
1882 {
1883  int el = 0;
1884  int eg = 0;
1885  int ind[8];
1886  NURBSPatchMap p2g(this);
1887  KnotVector *kv[3];
1888 
1889  for (int p = 0; p < GetNP(); p++)
1890  {
1891  p2g.SetPatchVertexMap(p, kv);
1892  int nx = p2g.nx();
1893  int ny = p2g.ny();
1894  int nz = p2g.nz();
1895 
1896  int patch_attr = patchTopo->GetAttribute(p);
1897 
1898  for (int k = 0; k < nz; k++)
1899  {
1900  for (int j = 0; j < ny; j++)
1901  {
1902  for (int i = 0; i < nx; i++)
1903  {
1904  if (activeElem[eg])
1905  {
1906  ind[0] = activeVert[p2g(i, j, k)];
1907  ind[1] = activeVert[p2g(i+1,j, k)];
1908  ind[2] = activeVert[p2g(i+1,j+1,k)];
1909  ind[3] = activeVert[p2g(i, j+1,k)];
1910 
1911  ind[4] = activeVert[p2g(i, j, k+1)];
1912  ind[5] = activeVert[p2g(i+1,j, k+1)];
1913  ind[6] = activeVert[p2g(i+1,j+1,k+1)];
1914  ind[7] = activeVert[p2g(i, j+1,k+1)];
1915 
1916  elements[el] = new Hexahedron(ind, patch_attr);
1917  el++;
1918  }
1919  eg++;
1920  }
1921  }
1922  }
1923  }
1924 }
1925 
1927 {
1928  boundary.SetSize(GetNBE());
1929 
1930  if (Dimension() == 2)
1931  {
1932  Get2DBdrElementTopo(boundary);
1933  }
1934  else
1935  {
1936  Get3DBdrElementTopo(boundary);
1937  }
1938 }
1939 
1941 {
1942  int g_be, l_be;
1943  int ind[2], okv[1];
1944  NURBSPatchMap p2g(this);
1945  KnotVector *kv[1];
1946 
1947  g_be = l_be = 0;
1948  for (int b = 0; b < GetNBP(); b++)
1949  {
1950  p2g.SetBdrPatchVertexMap(b, kv, okv);
1951  int nx = p2g.nx();
1952 
1953  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
1954 
1955  for (int i = 0; i < nx; i++)
1956  {
1957  if (activeBdrElem[g_be])
1958  {
1959  int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
1960  ind[0] = activeVert[p2g[_i ]];
1961  ind[1] = activeVert[p2g[_i+1]];
1962 
1963  boundary[l_be] = new Segment(ind, bdr_patch_attr);
1964  l_be++;
1965  }
1966  g_be++;
1967  }
1968  }
1969 }
1970 
1972 {
1973  int g_be, l_be;
1974  int ind[4], okv[2];
1975  NURBSPatchMap p2g(this);
1976  KnotVector *kv[2];
1977 
1978  g_be = l_be = 0;
1979  for (int b = 0; b < GetNBP(); b++)
1980  {
1981  p2g.SetBdrPatchVertexMap(b, kv, okv);
1982  int nx = p2g.nx();
1983  int ny = p2g.ny();
1984 
1985  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
1986 
1987  for (int j = 0; j < ny; j++)
1988  {
1989  int _j = (okv[1] >= 0) ? j : (ny - 1 - j);
1990  for (int i = 0; i < nx; i++)
1991  {
1992  if (activeBdrElem[g_be])
1993  {
1994  int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
1995  ind[0] = activeVert[p2g(_i, _j )];
1996  ind[1] = activeVert[p2g(_i+1,_j )];
1997  ind[2] = activeVert[p2g(_i+1,_j+1)];
1998  ind[3] = activeVert[p2g(_i, _j+1)];
1999 
2000  boundary[l_be] = new Quadrilateral(ind, bdr_patch_attr);
2001  l_be++;
2002  }
2003  g_be++;
2004  }
2005  }
2006  }
2007 }
2008 
2010 {
2011  activeDof.SetSize(GetNTotalDof());
2012  activeDof = 0;
2013 
2014  if (Dimension() == 2)
2015  {
2016  Generate2DElementDofTable();
2017  }
2018  else
2019  {
2020  Generate3DElementDofTable();
2021  }
2022 
2023  NumOfActiveDofs = 0;
2024  for (int d = 0; d < GetNTotalDof(); d++)
2025  if (activeDof[d])
2026  {
2027  NumOfActiveDofs++;
2028  activeDof[d] = NumOfActiveDofs;
2029  }
2030 
2031  int *dof = el_dof->GetJ();
2032  int ndof = el_dof->Size_of_connections();
2033  for (int i = 0; i < ndof; i++)
2034  {
2035  dof[i] = activeDof[dof[i]] - 1;
2036  }
2037 }
2038 
2040 {
2041  int el = 0;
2042  int eg = 0;
2043  KnotVector *kv[2];
2044  NURBSPatchMap p2g(this);
2045 
2046  el_dof = new Table(NumOfActiveElems, (Order + 1)*(Order + 1));
2047  el_to_patch.SetSize(NumOfActiveElems);
2048  el_to_IJK.SetSize(NumOfActiveElems, 2);
2049 
2050  for (int p = 0; p < GetNP(); p++)
2051  {
2052  p2g.SetPatchDofMap(p, kv);
2053 
2054  // Load dofs
2055  for (int j = 0; j < kv[1]->GetNKS(); j++)
2056  {
2057  if (kv[1]->isElement(j))
2058  {
2059  for (int i = 0; i < kv[0]->GetNKS(); i++)
2060  {
2061  if (kv[0]->isElement(i))
2062  {
2063  if (activeElem[eg])
2064  {
2065  int *dofs = el_dof->GetRow(el);
2066  int idx = 0;
2067  for (int jj = 0; jj <= Order; jj++)
2068  {
2069  for (int ii = 0; ii <= Order; ii++)
2070  {
2071  dofs[idx] = p2g(i+ii,j+jj);
2072  activeDof[dofs[idx]] = 1;
2073  idx++;
2074  }
2075  }
2076  el_to_patch[el] = p;
2077  el_to_IJK(el,0) = i;
2078  el_to_IJK(el,1) = j;
2079 
2080  el++;
2081  }
2082  eg++;
2083  }
2084  }
2085  }
2086  }
2087  }
2088 }
2089 
2091 {
2092  int el = 0;
2093  int eg = 0;
2094  int idx;
2095  KnotVector *kv[3];
2096  NURBSPatchMap p2g(this);
2097 
2098  el_dof = new Table(NumOfActiveElems, (Order + 1)*(Order + 1)*(Order + 1));
2099  el_to_patch.SetSize(NumOfActiveElems);
2100  el_to_IJK.SetSize(NumOfActiveElems, 3);
2101 
2102  for (int p = 0; p < GetNP(); p++)
2103  {
2104  p2g.SetPatchDofMap(p, kv);
2105 
2106  // Load dofs
2107  for (int k = 0; k < kv[2]->GetNKS(); k++)
2108  {
2109  if (kv[2]->isElement(k))
2110  {
2111  for (int j = 0; j < kv[1]->GetNKS(); j++)
2112  {
2113  if (kv[1]->isElement(j))
2114  {
2115  for (int i = 0; i < kv[0]->GetNKS(); i++)
2116  {
2117  if (kv[0]->isElement(i))
2118  {
2119  if (activeElem[eg])
2120  {
2121  idx = 0;
2122  int *dofs = el_dof->GetRow(el);
2123  for (int kk = 0; kk <= Order; kk++)
2124  {
2125  for (int jj = 0; jj <= Order; jj++)
2126  {
2127  for (int ii = 0; ii <= Order; ii++)
2128  {
2129  dofs[idx] = p2g(i+ii,j+jj,k+kk);
2130  activeDof[dofs[idx]] = 1;
2131  idx++;
2132  }
2133  }
2134  }
2135 
2136  el_to_patch[el] = p;
2137  el_to_IJK(el,0) = i;
2138  el_to_IJK(el,1) = j;
2139  el_to_IJK(el,2) = k;
2140 
2141  el++;
2142  }
2143  eg++;
2144  }
2145  }
2146  }
2147  }
2148  }
2149  }
2150  }
2151 }
2152 
2154 {
2155  if (Dimension() == 2)
2156  {
2157  Generate2DBdrElementDofTable();
2158  }
2159  else
2160  {
2161  Generate3DBdrElementDofTable();
2162  }
2163 
2164  int *dof = bel_dof->GetJ();
2165  int ndof = bel_dof->Size_of_connections();
2166  for (int i = 0; i < ndof; i++)
2167  {
2168  dof[i] = activeDof[dof[i]] - 1;
2169  }
2170 }
2171 
2173 {
2174  int gbe = 0;
2175  int lbe = 0, okv[1];
2176  KnotVector *kv[1];
2177  NURBSPatchMap p2g(this);
2178 
2179  bel_dof = new Table(NumOfActiveBdrElems, Order + 1);
2180  bel_to_patch.SetSize(NumOfActiveBdrElems);
2181  bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
2182 
2183  for (int b = 0; b < GetNBP(); b++)
2184  {
2185  p2g.SetBdrPatchDofMap(b, kv, okv);
2186  int nx = p2g.nx(); // NCP-1
2187 
2188  // Load dofs
2189  int nks0 = kv[0]->GetNKS();
2190  for (int i = 0; i < nks0; i++)
2191  {
2192  if (kv[0]->isElement(i))
2193  {
2194  if (activeBdrElem[gbe])
2195  {
2196  int *dofs = bel_dof->GetRow(lbe);
2197  int idx = 0;
2198  for (int ii = 0; ii <= Order; ii++)
2199  {
2200  dofs[idx] = p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)];
2201  idx++;
2202  }
2203  bel_to_patch[lbe] = b;
2204  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2205  lbe++;
2206  }
2207  gbe++;
2208  }
2209  }
2210  }
2211 }
2212 
2214 {
2215  int gbe = 0;
2216  int lbe = 0, okv[2];
2217  KnotVector *kv[2];
2218  NURBSPatchMap p2g(this);
2219 
2220  bel_dof = new Table(NumOfActiveBdrElems, (Order + 1)*(Order + 1));
2221  bel_to_patch.SetSize(NumOfActiveBdrElems);
2222  bel_to_IJK.SetSize(NumOfActiveBdrElems, 2);
2223 
2224  for (int b = 0; b < GetNBP(); b++)
2225  {
2226  p2g.SetBdrPatchDofMap(b, kv, okv);
2227  int nx = p2g.nx(); // NCP0-1
2228  int ny = p2g.ny(); // NCP1-1
2229 
2230  // Load dofs
2231  int nks0 = kv[0]->GetNKS();
2232  int nks1 = kv[1]->GetNKS();
2233  for (int j = 0; j < nks1; j++)
2234  {
2235  if (kv[1]->isElement(j))
2236  {
2237  for (int i = 0; i < nks0; i++)
2238  {
2239  if (kv[0]->isElement(i))
2240  {
2241  if (activeBdrElem[gbe])
2242  {
2243  int *dofs = bel_dof->GetRow(lbe);
2244  int idx = 0;
2245  for (int jj = 0; jj <= Order; jj++)
2246  {
2247  int _jj = (okv[1] >= 0) ? (j+jj) : (ny-j-jj);
2248  for (int ii = 0; ii <= Order; ii++)
2249  {
2250  int _ii = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
2251  dofs[idx] = p2g(_ii,_jj);
2252  idx++;
2253  }
2254  }
2255  bel_to_patch[lbe] = b;
2256  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2257  bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
2258  lbe++;
2259  }
2260  gbe++;
2261  }
2262  }
2263  }
2264  }
2265  }
2266 }
2267 
2269 {
2270  lvert_vert.SetSize(GetNV());
2271  for (int gv = 0; gv < GetGNV(); gv++)
2272  if (activeVert[gv] >= 0)
2273  {
2274  lvert_vert[activeVert[gv]] = gv;
2275  }
2276 }
2277 
2279 {
2280  lelem_elem.SetSize(GetNE());
2281  for (int le = 0, ge = 0; ge < GetGNE(); ge++)
2282  if (activeElem[ge])
2283  {
2284  lelem_elem[le++] = ge;
2285  }
2286 }
2287 
2289 {
2290  const NURBSFiniteElement *NURBSFE =
2291  dynamic_cast<const NURBSFiniteElement *>(FE);
2292 
2293  if (NURBSFE->GetElement() != i)
2294  {
2295  Array<int> dofs;
2296  NURBSFE->SetIJK(el_to_IJK.GetRow(i));
2297  if (el_to_patch[i] != NURBSFE->GetPatch())
2298  {
2299  GetPatchKnotVectors(el_to_patch[i], NURBSFE->KnotVectors());
2300  NURBSFE->SetPatch(el_to_patch[i]);
2301  }
2302  el_dof->GetRow(i, dofs);
2303  weights.GetSubVector(dofs, NURBSFE->Weights());
2304  NURBSFE->SetElement(i);
2305  }
2306 }
2307 
2309 {
2310  const NURBSFiniteElement *NURBSFE =
2311  dynamic_cast<const NURBSFiniteElement *>(BE);
2312 
2313  if (NURBSFE->GetElement() != i)
2314  {
2315  Array<int> dofs;
2316  NURBSFE->SetIJK(bel_to_IJK.GetRow(i));
2317  if (bel_to_patch[i] != NURBSFE->GetPatch())
2318  {
2319  GetBdrPatchKnotVectors(bel_to_patch[i], NURBSFE->KnotVectors());
2320  NURBSFE->SetPatch(bel_to_patch[i]);
2321  }
2322  bel_dof->GetRow(i, dofs);
2323  weights.GetSubVector(dofs, NURBSFE->Weights());
2324  NURBSFE->SetElement(i);
2325  }
2326 }
2327 
2329 {
2330  delete el_dof;
2331  delete bel_dof;
2332 
2333  if (patches.Size() == 0)
2334  {
2335  GetPatchNets(Nodes);
2336  }
2337 }
2338 
2340 {
2341  if (patches.Size() == 0) { return; }
2342 
2343  SetSolutionVector(Nodes);
2344  patches.SetSize(0);
2345 }
2346 
2348 {
2349  if (patches.Size() == 0)
2350  mfem_error("NURBSExtension::SetKnotsFromPatches :"
2351  " No patches available!");
2352 
2354 
2355  for (int p = 0; p < patches.Size(); p++)
2356  {
2357  GetPatchKnotVectors(p, kv);
2358 
2359  for (int i = 0; i < kv.Size(); i++)
2360  {
2361  *kv[i] = *patches[p]->GetKV(i);
2362  }
2363  }
2364 
2365  Order = knotVectors[0]->GetOrder();
2366  for (int i = 0; i < NumOfKnotVectors; i++)
2367  {
2368  if (Order != knotVectors[i]->GetOrder())
2369  mfem_error("NURBSExtension::Reset :\n"
2370  " Variable orders are not supported!");
2371  }
2372 
2373  GenerateOffsets();
2374  CountElements();
2375  CountBdrElements();
2376 
2377  // all elements must be active
2378  NumOfActiveElems = NumOfElements;
2379  activeElem.SetSize(NumOfElements);
2380  activeElem = true;
2381 
2382  GenerateActiveVertices();
2383  GenerateElementDofTable();
2384  GenerateActiveBdrElems();
2385  GenerateBdrElementDofTable();
2386 }
2387 
2389 {
2390  for (int p = 0; p < patches.Size(); p++)
2391  {
2392  patches[p]->DegreeElevate(t);
2393  }
2394 }
2395 
2397 {
2398  for (int p = 0; p < patches.Size(); p++)
2399  {
2400  patches[p]->UniformRefinement();
2401  }
2402 }
2403 
2405 {
2406  Array<int> edges;
2407  Array<int> orient;
2408 
2409  Array<KnotVector *> pkv(Dimension());
2410 
2411  for (int p = 0; p < patches.Size(); p++)
2412  {
2413  patchTopo->GetElementEdges(p, edges, orient);
2414 
2415  if (Dimension()==2)
2416  {
2417  pkv[0] = kv[KnotInd(edges[0])];
2418  pkv[1] = kv[KnotInd(edges[1])];
2419  }
2420  else
2421  {
2422  pkv[0] = kv[KnotInd(edges[0])];
2423  pkv[1] = kv[KnotInd(edges[3])];
2424  pkv[2] = kv[KnotInd(edges[8])];
2425  }
2426 
2427  patches[p]->KnotInsert(pkv);
2428  }
2429 }
2430 
2432 {
2433  if (Dimension() == 2)
2434  {
2435  Get2DPatchNets(coords);
2436  }
2437  else
2438  {
2439  Get3DPatchNets(coords);
2440  }
2441 }
2442 
2444 {
2445  Array<KnotVector *> kv(2);
2446  NURBSPatchMap p2g(this);
2447 
2448  patches.SetSize(GetNP());
2449  for (int p = 0; p < GetNP(); p++)
2450  {
2451  p2g.SetPatchDofMap(p, kv);
2452  patches[p] = new NURBSPatch(kv, 2+1);
2453  NURBSPatch &Patch = *patches[p];
2454 
2455  for (int j = 0; j < kv[1]->GetNCP(); j++)
2456  {
2457  for (int i = 0; i < kv[0]->GetNCP(); i++)
2458  {
2459  int l = p2g(i,j);
2460  for (int d = 0; d < 2; d++)
2461  {
2462  Patch(i,j,d) = coords(l*2 + d)*weights(l);
2463  }
2464  Patch(i,j,2) = weights(l);
2465  }
2466  }
2467  }
2468 }
2469 
2471 {
2472  Array<KnotVector *> kv(3);
2473  NURBSPatchMap p2g(this);
2474 
2475  patches.SetSize(GetNP());
2476  for (int p = 0; p < GetNP(); p++)
2477  {
2478  p2g.SetPatchDofMap(p, kv);
2479  patches[p] = new NURBSPatch(kv, 3+1);
2480  NURBSPatch &Patch = *patches[p];
2481 
2482  for (int k = 0; k < kv[2]->GetNCP(); k++)
2483  {
2484  for (int j = 0; j < kv[1]->GetNCP(); j++)
2485  {
2486  for (int i = 0; i < kv[0]->GetNCP(); i++)
2487  {
2488  int l = p2g(i,j,k);
2489  for (int d = 0; d < 3; d++)
2490  {
2491  Patch(i,j,k,d) = coords(l*3 + d)*weights(l);
2492  }
2493  Patch(i,j,k,3) = weights(l);
2494  }
2495  }
2496  }
2497  }
2498 }
2499 
2501 {
2502  if (Dimension() == 2)
2503  {
2504  Set2DSolutionVector(coords);
2505  }
2506  else
2507  {
2508  Set3DSolutionVector(coords);
2509  }
2510 }
2511 
2513 {
2514  Array<KnotVector *> kv(2);
2515  NURBSPatchMap p2g(this);
2516 
2517  weights.SetSize(GetNDof());
2518  for (int p = 0; p < GetNP(); p++)
2519  {
2520  p2g.SetPatchDofMap(p, kv);
2521  NURBSPatch &Patch = *patches[p];
2522 
2523  for (int j = 0; j < kv[1]->GetNCP(); j++)
2524  {
2525  for (int i = 0; i < kv[0]->GetNCP(); i++)
2526  {
2527  int l = p2g(i,j);
2528  for (int d = 0; d < 2; d++)
2529  {
2530  coords(l*2 + d) = Patch(i,j,d)/Patch(i,j,2);
2531  }
2532  weights(l) = Patch(i,j,2);
2533  }
2534  }
2535  delete patches[p];
2536  }
2537 }
2538 
2540 {
2541  Array<KnotVector *> kv(3);
2542  NURBSPatchMap p2g(this);
2543 
2544  weights.SetSize(GetNDof());
2545  for (int p = 0; p < GetNP(); p++)
2546  {
2547  p2g.SetPatchDofMap(p, kv);
2548  NURBSPatch &Patch = *patches[p];
2549 
2550  for (int k = 0; k < kv[2]->GetNCP(); k++)
2551  {
2552  for (int j = 0; j < kv[1]->GetNCP(); j++)
2553  {
2554  for (int i = 0; i < kv[0]->GetNCP(); i++)
2555  {
2556  int l = p2g(i,j,k);
2557  for (int d = 0; d < 3; d++)
2558  {
2559  coords(l*3 + d) = Patch(i,j,k,d)/Patch(i,j,k,3);
2560  }
2561  weights(l) = Patch(i,j,k,3);
2562  }
2563  }
2564  }
2565  delete patches[p];
2566  }
2567 }
2568 
2569 
2570 #ifdef MFEM_USE_MPI
2572  int *part, const Array<bool> &active_bel)
2573  : gtopo(comm)
2574 {
2575  if (parent->NumOfActiveElems < parent->NumOfElements)
2576  // SetActive (BuildGroups?) and the way the weights are copied
2577  // do not support this case
2578  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
2579  " all elements in the parent must be active!");
2580 
2581  patchTopo = parent->patchTopo;
2582  // steal ownership of patchTopo from the 'parent' NURBS extension
2583  if (!parent->own_topo)
2584  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
2585  " parent does not own the patch topology!");
2586  own_topo = 1;
2587  parent->own_topo = 0;
2588 
2589  parent->edge_to_knot.Copy(edge_to_knot);
2590 
2591  Order = parent->GetOrder();
2592 
2593  NumOfKnotVectors = parent->GetNKV();
2594  knotVectors.SetSize(NumOfKnotVectors);
2595  for (int i = 0; i < NumOfKnotVectors; i++)
2596  {
2597  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
2598  }
2599 
2600  GenerateOffsets();
2601  CountElements();
2602  CountBdrElements();
2603 
2604  // copy 'part' to 'partitioning'
2605  partitioning = new int[GetGNE()];
2606  for (int i = 0; i < GetGNE(); i++)
2607  {
2608  partitioning[i] = part[i];
2609  }
2610  SetActive(partitioning, active_bel);
2611 
2614  // GenerateActiveBdrElems(); // done by SetActive for now
2616 
2617  Table *serial_elem_dof = parent->GetElementDofTable();
2618  BuildGroups(partitioning, *serial_elem_dof);
2619 
2620  weights.SetSize(GetNDof());
2621  // copy weights from parent
2622  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
2623  {
2624  if (activeElem[gel])
2625  {
2626  int ndofs = el_dof->RowSize(lel);
2627  int *ldofs = el_dof->GetRow(lel);
2628  int *gdofs = serial_elem_dof->GetRow(gel);
2629  for (int i = 0; i < ndofs; i++)
2630  {
2631  weights(ldofs[i]) = parent->weights(gdofs[i]);
2632  }
2633  lel++;
2634  }
2635  }
2636 }
2637 
2639  ParNURBSExtension *par_parent)
2640  : gtopo(par_parent->gtopo.GetComm())
2641 {
2642  // steal all data from parent
2643  Order = parent->Order;
2644 
2645  patchTopo = parent->patchTopo;
2646  own_topo = parent->own_topo;
2647  parent->own_topo = 0;
2648 
2649  Swap(edge_to_knot, parent->edge_to_knot);
2650 
2652  Swap(knotVectors, parent->knotVectors);
2653 
2654  NumOfVertices = parent->NumOfVertices;
2655  NumOfElements = parent->NumOfElements;
2657  NumOfDofs = parent->NumOfDofs;
2658 
2659  Swap(v_meshOffsets, parent->v_meshOffsets);
2660  Swap(e_meshOffsets, parent->e_meshOffsets);
2661  Swap(f_meshOffsets, parent->f_meshOffsets);
2662  Swap(p_meshOffsets, parent->p_meshOffsets);
2663 
2668 
2672  NumOfActiveDofs = parent->NumOfActiveDofs;
2673 
2674  Swap(activeVert, parent->activeVert);
2675  Swap(activeElem, parent->activeElem);
2676  Swap(activeBdrElem, parent->activeBdrElem);
2677  Swap(activeDof, parent->activeDof);
2678 
2679  el_dof = parent->el_dof;
2680  bel_dof = parent->bel_dof;
2681  parent->el_dof = parent->bel_dof = NULL;
2682 
2683  Swap(el_to_patch, parent->el_to_patch);
2684  Swap(bel_to_patch, parent->bel_to_patch);
2685  Swap(el_to_IJK, parent->el_to_IJK);
2686  Swap(bel_to_IJK, parent->bel_to_IJK);
2687 
2688  Swap(weights, parent->weights);
2689 
2690  delete parent;
2691 
2692  partitioning = NULL;
2693 
2694 #ifdef MFEM_DEBUG
2695  if (par_parent->partitioning == NULL)
2696  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
2697  " parent ParNURBSExtension has no partitioning!");
2698 #endif
2699 
2700  Table *serial_elem_dof = GetGlobalElementDofTable();
2701  BuildGroups(par_parent->partitioning, *serial_elem_dof);
2702  delete serial_elem_dof;
2703 }
2704 
2705 Table *ParNURBSExtension::GetGlobalElementDofTable()
2706 {
2707  if (Dimension() == 2)
2708  {
2709  return Get2DGlobalElementDofTable();
2710  }
2711  else
2712  {
2713  return Get3DGlobalElementDofTable();
2714  }
2715 }
2716 
2717 Table *ParNURBSExtension::Get2DGlobalElementDofTable()
2718 {
2719  int el = 0;
2720  KnotVector *kv[2];
2721  NURBSPatchMap p2g(this);
2722 
2723  Table *gel_dof = new Table(GetGNE(), (Order + 1)*(Order + 1));
2724 
2725  for (int p = 0; p < GetNP(); p++)
2726  {
2727  p2g.SetPatchDofMap(p, kv);
2728 
2729  // Load dofs
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  int *dofs = gel_dof->GetRow(el);
2739  int idx = 0;
2740  for (int jj = 0; jj <= Order; jj++)
2741  {
2742  for (int ii = 0; ii <= Order; ii++)
2743  {
2744  dofs[idx] = p2g(i+ii,j+jj);
2745  idx++;
2746  }
2747  }
2748  el++;
2749  }
2750  }
2751  }
2752  }
2753  }
2754  return gel_dof;
2755 }
2756 
2757 Table *ParNURBSExtension::Get3DGlobalElementDofTable()
2758 {
2759  int el = 0;
2760  KnotVector *kv[3];
2761  NURBSPatchMap p2g(this);
2762 
2763  Table *gel_dof = new Table(GetGNE(), (Order + 1)*(Order + 1)*(Order + 1));
2764 
2765  for (int p = 0; p < GetNP(); p++)
2766  {
2767  p2g.SetPatchDofMap(p, kv);
2768 
2769  // Load dofs
2770  for (int k = 0; k < kv[2]->GetNKS(); k++)
2771  {
2772  if (kv[2]->isElement(k))
2773  {
2774  for (int j = 0; j < kv[1]->GetNKS(); j++)
2775  {
2776  if (kv[1]->isElement(j))
2777  {
2778  for (int i = 0; i < kv[0]->GetNKS(); i++)
2779  {
2780  if (kv[0]->isElement(i))
2781  {
2782  int *dofs = gel_dof->GetRow(el);
2783  int idx = 0;
2784  for (int kk = 0; kk <= Order; kk++)
2785  {
2786  for (int jj = 0; jj <= Order; jj++)
2787  {
2788  for (int ii = 0; ii <= Order; ii++)
2789  {
2790  dofs[idx] = p2g(i+ii,j+jj,k+kk);
2791  idx++;
2792  }
2793  }
2794  }
2795  el++;
2796  }
2797  }
2798  }
2799  }
2800  }
2801  }
2802  }
2803  return gel_dof;
2804 }
2805 
2806 void ParNURBSExtension::SetActive(int *_partitioning,
2807  const Array<bool> &active_bel)
2808 {
2810  activeElem = false;
2811  NumOfActiveElems = 0;
2812  int MyRank = gtopo.MyRank();
2813  for (int i = 0; i < GetGNE(); i++)
2814  if (_partitioning[i] == MyRank)
2815  {
2816  activeElem[i] = true;
2817  NumOfActiveElems++;
2818  }
2819 
2820  active_bel.Copy(activeBdrElem);
2821  NumOfActiveBdrElems = 0;
2822  for (int i = 0; i < GetGNBE(); i++)
2823  if (activeBdrElem[i])
2824  {
2826  }
2827 }
2828 
2829 void ParNURBSExtension::BuildGroups(int *_partitioning, const Table &elem_dof)
2830 {
2831  Table dof_proc;
2832 
2833  ListOfIntegerSets groups;
2834  IntegerSet group;
2835 
2836  Transpose(elem_dof, dof_proc); // dof_proc is dof_elem
2837  // convert elements to processors
2838  for (int i = 0; i < dof_proc.Size_of_connections(); i++)
2839  {
2840  dof_proc.GetJ()[i] = _partitioning[dof_proc.GetJ()[i]];
2841  }
2842 
2843  // the first group is the local one
2844  int MyRank = gtopo.MyRank();
2845  group.Recreate(1, &MyRank);
2846  groups.Insert(group);
2847 
2848  int dof = 0;
2850  for (int d = 0; d < GetNTotalDof(); d++)
2851  if (activeDof[d])
2852  {
2853  group.Recreate(dof_proc.RowSize(d), dof_proc.GetRow(d));
2854  ldof_group[dof] = groups.Insert(group);
2855 
2856  dof++;
2857  }
2858 
2859  gtopo.Create(groups, 1822);
2860 }
2861 #endif
2862 
2863 
2864 void NURBSPatchMap::GetPatchKnotVectors(int p, KnotVector *kv[])
2865 {
2866  Ext->patchTopo->GetElementVertices(p, verts);
2867  Ext->patchTopo->GetElementEdges(p, edges, oedge);
2868  if (Ext->Dimension() == 2)
2869  {
2870  kv[0] = Ext->KnotVec(edges[0]);
2871  kv[1] = Ext->KnotVec(edges[1]);
2872  }
2873  else
2874  {
2875  Ext->patchTopo->GetElementFaces(p, faces, oface);
2876 
2877  kv[0] = Ext->KnotVec(edges[0]);
2878  kv[1] = Ext->KnotVec(edges[3]);
2879  kv[2] = Ext->KnotVec(edges[8]);
2880  }
2881  opatch = 0;
2882 }
2883 
2884 void NURBSPatchMap::GetBdrPatchKnotVectors(int p, KnotVector *kv[], int *okv)
2885 {
2886  Ext->patchTopo->GetBdrElementVertices(p, verts);
2887  Ext->patchTopo->GetBdrElementEdges(p, edges, oedge);
2888  kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]);
2889  if (Ext->Dimension() == 3)
2890  {
2891  faces.SetSize(1);
2892  Ext->patchTopo->GetBdrElementFace(p, &faces[0], &opatch);
2893  kv[1] = Ext->KnotVec(edges[1], oedge[1], &okv[1]);
2894  }
2895  else
2896  {
2897  opatch = oedge[0];
2898  }
2899 }
2900 
2902 {
2903  GetPatchKnotVectors(p, kv);
2904 
2905  I = kv[0]->GetNE() - 1;
2906  J = kv[1]->GetNE() - 1;
2907 
2908  for (int i = 0; i < verts.Size(); i++)
2909  {
2910  verts[i] = Ext->v_meshOffsets[verts[i]];
2911  }
2912 
2913  for (int i = 0; i < edges.Size(); i++)
2914  {
2915  edges[i] = Ext->e_meshOffsets[edges[i]];
2916  }
2917 
2918  if (Ext->Dimension() == 3)
2919  {
2920  K = kv[2]->GetNE() - 1;
2921 
2922  for (int i = 0; i < faces.Size(); i++)
2923  {
2924  faces[i] = Ext->f_meshOffsets[faces[i]];
2925  }
2926  }
2927 
2928  pOffset = Ext->p_meshOffsets[p];
2929 }
2930 
2932 {
2933  GetPatchKnotVectors(p, kv);
2934 
2935  I = kv[0]->GetNCP() - 2;
2936  J = kv[1]->GetNCP() - 2;
2937 
2938  for (int i = 0; i < verts.Size(); i++)
2939  {
2940  verts[i] = Ext->v_spaceOffsets[verts[i]];
2941  }
2942 
2943  for (int i = 0; i < edges.Size(); i++)
2944  {
2945  edges[i] = Ext->e_spaceOffsets[edges[i]];
2946  }
2947 
2948  if (Ext->Dimension() == 3)
2949  {
2950  K = kv[2]->GetNCP() - 2;
2951 
2952  for (int i = 0; i < faces.Size(); i++)
2953  {
2954  faces[i] = Ext->f_spaceOffsets[faces[i]];
2955  }
2956  }
2957 
2958  pOffset = Ext->p_spaceOffsets[p];
2959 }
2960 
2962 {
2963  GetBdrPatchKnotVectors(p, kv, okv);
2964 
2965  I = kv[0]->GetNE() - 1;
2966 
2967  for (int i = 0; i < verts.Size(); i++)
2968  {
2969  verts[i] = Ext->v_meshOffsets[verts[i]];
2970  }
2971 
2972  if (Ext->Dimension() == 2)
2973  {
2974  pOffset = Ext->e_meshOffsets[edges[0]];
2975  }
2976  else
2977  {
2978  J = kv[1]->GetNE() - 1;
2979 
2980  for (int i = 0; i < edges.Size(); i++)
2981  {
2982  edges[i] = Ext->e_meshOffsets[edges[i]];
2983  }
2984 
2985  pOffset = Ext->f_meshOffsets[faces[0]];
2986  }
2987 }
2988 
2989 void NURBSPatchMap::SetBdrPatchDofMap(int p, KnotVector *kv[], int *okv)
2990 {
2991  GetBdrPatchKnotVectors(p, kv, okv);
2992 
2993  I = kv[0]->GetNCP() - 2;
2994 
2995  for (int i = 0; i < verts.Size(); i++)
2996  {
2997  verts[i] = Ext->v_spaceOffsets[verts[i]];
2998  }
2999 
3000  if (Ext->Dimension() == 2)
3001  {
3002  pOffset = Ext->e_spaceOffsets[edges[0]];
3003  }
3004  else
3005  {
3006  J = kv[1]->GetNCP() - 2;
3007 
3008  for (int i = 0; i < edges.Size(); i++)
3009  {
3010  edges[i] = Ext->e_spaceOffsets[edges[i]];
3011  }
3012 
3013  pOffset = Ext->f_spaceOffsets[faces[0]];
3014  }
3015 }
3016 
3017 }
Array< KnotVector * > knotVectors
Definition: nurbs.hpp:174
Abstract class for Finite Elements.
Definition: fe.hpp:44
Array< KnotVector * > kv
Definition: nurbs.hpp:89
Array< int > f_meshOffsets
Definition: nurbs.hpp:180
KnotVector * DegreeElevate(int t) const
Definition: nurbs.cpp:54
void Get3DPatchNets(const Vector &Nodes)
Definition: nurbs.cpp:2470
void Set3DSolutionVector(Vector &Nodes)
Definition: nurbs.cpp:2539
void Create(ListOfIntegerSets &groups, int mpitag)
NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP)
Definition: nurbs.cpp:388
Array2D< int > bel_to_IJK
Definition: nurbs.hpp:194
int Size() const
Logical size of the array.
Definition: array.hpp:109
void SetCoordsFromPatches(Vector &Nodes)
Definition: nurbs.cpp:2339
int Size() const
Definition: nurbs.hpp:50
Array< int > activeVert
Definition: nurbs.hpp:166
KnotVector * KnotVec(int edge)
Definition: nurbs.hpp:478
void swap(NURBSPatch *np)
Definition: nurbs.cpp:403
void Get2DBdrElementTopo(Array< Element * > &boundary)
Definition: nurbs.cpp:1940
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
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:3372
void DegreeElevate(int dir, int t)
Definition: nurbs.cpp:681
int NumOfElements
Definition: nurbs.hpp:35
NURBSPatch * Interpolate(NURBSPatch &p1, NURBSPatch &p2)
Definition: nurbs.cpp:1052
void SetPatchDofMap(int p, KnotVector *kv[])
Definition: nurbs.cpp:2931
void Generate2DElementDofTable()
Definition: nurbs.cpp:2039
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Definition: nurbs.cpp:1565
Array< int > ldof_group
Definition: nurbs.hpp:350
int GetPatch() const
Definition: fe.hpp:2075
void Generate3DElementDofTable()
Definition: nurbs.cpp:2090
void KnotInsert(int dir, const KnotVector &knot)
Definition: nurbs.cpp:543
void MergeWeights(Mesh *mesh_array[], int num_pieces)
Definition: nurbs.cpp:1540
void SetSize(int s)
Resize the vector if the new size is different.
Definition: vector.hpp:263
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:132
void UniformRefinement()
Definition: nurbs.cpp:2396
void GetElementFaces(int i, Array< int > &, Array< int > &) const
Return the indices and the orientations of all faces of element i.
Definition: mesh.cpp:3565
int GetNE() const
Definition: nurbs.hpp:46
void Print(std::ostream &out)
Definition: nurbs.cpp:442
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:160
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:457
void Rotate3D(double normal[], double angle)
Definition: nurbs.cpp:1007
Array< int > e_meshOffsets
Definition: nurbs.hpp:179
int GetNKS() const
Definition: nurbs.hpp:47
void GetElements()
Count the number of elements.
Definition: nurbs.cpp:95
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
int Size() const
Returns the size of the vector.
Definition: vector.hpp:86
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:179
void GetPatchNets(const Vector &Nodes)
Definition: nurbs.cpp:2431
void Generate3DBdrElementDofTable()
Definition: nurbs.cpp:2213
int NumOfControlPoints
Definition: nurbs.hpp:35
Array< int > e_spaceOffsets
Definition: nurbs.hpp:185
void GenerateBdrElementDofTable()
Definition: nurbs.cpp:2153
void GetElementVertices(int i, Array< int > &dofs) const
Returns the indices of the dofs of element i.
Definition: mesh.hpp:558
const KnotVector * GetKnotVector(int i) const
Definition: nurbs.hpp:300
void Get2DElementTopo(Array< Element * > &elements)
Definition: nurbs.cpp:1845
void GetElementLocalToGlobal(Array< int > &lelem_elem)
Definition: nurbs.cpp:2278
void init(int dim_)
Definition: nurbs.cpp:277
Data type quadrilateral element.
Array< int > edge_to_knot
Definition: nurbs.hpp:173
void GenerateActiveVertices()
Definition: nurbs.cpp:1451
Array< bool > activeElem
Definition: nurbs.hpp:167
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
Definition: nurbs.cpp:2268
Vector knot
Definition: nurbs.hpp:34
int GetElement() const
Definition: fe.hpp:2077
Array< int > el_to_patch
Definition: nurbs.hpp:191
void SetBdrPatchVertexMap(int p, KnotVector *kv[], int *okv)
Definition: nurbs.cpp:2961
void SwapDirections(int dir1, int dir2)
Definition: nurbs.cpp:944
friend class NURBSPatchMap
Definition: nurbs.hpp:155
void Get3DElementTopo(Array< Element * > &elements)
Definition: nurbs.cpp:1881
KnotVector()
Create KnotVector.
Definition: nurbs.hpp:39
void Get3DBdrElementTopo(Array< Element * > &boundary)
Definition: nurbs.cpp:1971
static const int MaxOrder
Definition: nurbs.hpp:32
Array< int > v_spaceOffsets
Definition: nurbs.hpp:184
void UniformRefinement(Vector &newknots) const
Definition: nurbs.cpp:81
void LoadFE(int i, const FiniteElement *FE)
Definition: nurbs.cpp:2288
void GetBdrElementTopo(Array< Element * > &boundary)
Definition: nurbs.cpp:1926
void GetElementTopo(Array< Element * > &elements)
Definition: nurbs.cpp:1831
void Generate2DBdrElementDofTable()
Definition: nurbs.cpp:2172
Data type hexahedron element.
Definition: hexahedron.hpp:22
int dim
Definition: ex3.cpp:47
void SetPatch(int p) const
Definition: fe.hpp:2076
void GetPatchKnotVectors(int p, Array< KnotVector * > &kv)
Definition: nurbs.cpp:1664
void GenerateActiveBdrElems()
Definition: nurbs.cpp:1518
double * data
Definition: nurbs.hpp:87
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:136
void LoadPatchTopo(std::istream &input, Array< int > &edge_to_knot)
Read NURBS patch/macro-element mesh.
Definition: mesh.cpp:2864
virtual ~NURBSExtension()
Destroy a NURBSExtension.
Definition: nurbs.cpp:1367
void GenerateOffsets()
Definition: nurbs.cpp:1704
Array< int > activeDof
Definition: nurbs.hpp:169
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:391
static void Get3DRotationMatrix(double n[], double angle, double r, DenseMatrix &T)
Definition: nurbs.cpp:968
void SetKnotsFromPatches()
Definition: nurbs.cpp:2347
void CalcShape(Vector &shape, int i, double xi)
Definition: nurbs.cpp:127
void UniformRefinement()
Definition: nurbs.cpp:525
Array< int > f_spaceOffsets
Definition: nurbs.hpp:186
Array< bool > activeBdrElem
Definition: nurbs.hpp:168
void SetData(double *d)
Definition: vector.hpp:69
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:263
void SetElement(int e) const
Definition: fe.hpp:2078
int MakeUniformDegree()
Definition: nurbs.cpp:1033
void GetBdrElementFace(int i, int *, int *) const
Return the index and the orientation of the face of bdr element i. (3D)
Definition: mesh.cpp:3603
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:322
Array< int > v_meshOffsets
Definition: nurbs.hpp:178
void mfem_error(const char *msg)
Definition: error.cpp:23
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:331
void GetBdrPatchKnotVectors(int p, Array< KnotVector * > &kv)
Definition: nurbs.cpp:1685
void SetSubVector(const Array< int > &dofs, const Vector &elemvect)
Definition: vector.cpp:493
Array< int > p_meshOffsets
Definition: nurbs.hpp:181
KnotVector & operator=(const KnotVector &kv)
Definition: nurbs.cpp:42
void Set2DSolutionVector(Vector &Nodes)
Definition: nurbs.cpp:2512
Table * GetElementDofTable()
Definition: nurbs.hpp:308
Array2D< int > el_to_IJK
Definition: nurbs.hpp:193
void CountBdrElements()
Definition: nurbs.cpp:1811
NURBSExtension * NURBSext
Definition: mesh.hpp:143
ParNURBSExtension(MPI_Comm comm, NURBSExtension *parent, int *partitioning, const Array< bool > &active_bel)
Definition: nurbs.cpp:2571
Array< KnotVector * > & KnotVectors() const
Definition: fe.hpp:2079
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:3350
NURBSPatch * Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
Definition: nurbs.cpp:1099
Array< int > p_spaceOffsets
Definition: nurbs.hpp:187
int findKnotSpan(double u) const
Definition: nurbs.cpp:214
int GetNCP() const
Definition: nurbs.hpp:48
void GenerateElementDofTable()
Definition: nurbs.cpp:2009
static void skip_comment_lines(std::istream &is, const char comment_char)
Definition: mesh.hpp:177
void GetBdrElementVertices(int i, Array< int > &dofs) const
Returns the indices of the dofs of boundary element i.
Definition: mesh.hpp:562
void SetIJK(int *IJK) const
Definition: fe.hpp:2074
void SetBdrPatchDofMap(int p, KnotVector *kv[], int *okv)
Definition: nurbs.cpp:2989
void Difference(const KnotVector &kv, Vector &diff) const
Definition: nurbs.cpp:243
Vector & Weights() const
Definition: fe.hpp:2080
void CalcDShape(Vector &grad, int i, double xi)
Definition: nurbs.cpp:158
void DegreeElevate(int t)
Definition: nurbs.cpp:2388
void Get2DPatchNets(const Vector &Nodes)
Definition: nurbs.cpp:2443
Vector data type.
Definition: vector.hpp:33
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:142
void LoadBE(int i, const FiniteElement *BE)
Definition: nurbs.cpp:2308
void Print(std::ostream &out) const
Definition: nurbs.cpp:120
int RowSize(int i) const
Definition: table.hpp:102
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:71
void PrintCharacteristics(std::ostream &out)
Definition: nurbs.cpp:1426
void CheckBdrPatches()
Definition: nurbs.cpp:1636
void SetPatchVertexMap(int p, KnotVector *kv[])
Definition: nurbs.cpp:2901
void FlipDirection(int dir)
Definition: nurbs.cpp:932
void SetSolutionVector(Vector &Nodes)
Definition: nurbs.cpp:2500
GroupTopology gtopo
Definition: nurbs.hpp:348
Data type line segment element.
Definition: segment.hpp:22
int GetOrder() const
Definition: nurbs.hpp:49
Array< int > bel_to_patch
Definition: nurbs.hpp:192
void ConvertToPatches(const Vector &Nodes)
Definition: nurbs.cpp:2328
void KnotInsert(Array< KnotVector * > &kv)
Definition: nurbs.cpp:2404
void Print(std::ostream &out) const
Definition: nurbs.cpp:1391
int SetLoopDirection(int dir)
Definition: nurbs.cpp:466