MFEM  v3.1 Finite element discretization library
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
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  }
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
1160 // from mesh.cpp
1161 extern void skip_comment_lines(std::istream &, const char);
1162
1164 {
1166  patchTopo = new Mesh;
1168  own_topo = 1;
1169
1170  CheckPatches();
1171  // CheckBdrPatches();
1172
1173  skip_comment_lines(input, '#');
1174
1175  // Read knotvectors or patches
1176  string ident;
1177  input >> ws >> ident; // 'knotvectors' or 'patches'
1178  if (ident == "knotvectors")
1179  {
1180  input >> NumOfKnotVectors;
1181  knotVectors.SetSize(NumOfKnotVectors);
1182  for (int i = 0; i < NumOfKnotVectors; i++)
1183  {
1184  knotVectors[i] = new KnotVector(input);
1185  if (knotVectors[i]->GetOrder() != knotVectors[0]->GetOrder())
1186  mfem_error("NURBSExtension::NURBSExtension :\n"
1187  " Variable orders are not supported!");
1188  }
1189  Order = knotVectors[0]->GetOrder();
1190  }
1191  else if (ident == "patches")
1192  {
1193  patches.SetSize(GetNP());
1194  for (int p = 0; p < patches.Size(); p++)
1195  {
1196  skip_comment_lines(input, '#');
1197  patches[p] = new NURBSPatch(input);
1198  }
1199
1200  NumOfKnotVectors = 0;
1201  for (int i = 0; i < patchTopo->GetNEdges(); i++)
1202  if (NumOfKnotVectors < KnotInd(i))
1203  {
1204  NumOfKnotVectors = KnotInd(i);
1205  }
1206  NumOfKnotVectors++;
1207  knotVectors.SetSize(NumOfKnotVectors);
1208  knotVectors = NULL;
1209
1210  Array<int> edges, oedge;
1211  for (int p = 0; p < patches.Size(); p++)
1212  {
1213  patchTopo->GetElementEdges(p, edges, oedge);
1214  if (Dimension() == 2)
1215  {
1216  if (knotVectors[KnotInd(edges[0])] == NULL)
1217  knotVectors[KnotInd(edges[0])] =
1218  new KnotVector(*patches[p]->GetKV(0));
1219  if (knotVectors[KnotInd(edges[1])] == NULL)
1220  knotVectors[KnotInd(edges[1])] =
1221  new KnotVector(*patches[p]->GetKV(1));
1222  }
1223  else
1224  {
1225  if (knotVectors[KnotInd(edges[0])] == NULL)
1226  knotVectors[KnotInd(edges[0])] =
1227  new KnotVector(*patches[p]->GetKV(0));
1228  if (knotVectors[KnotInd(edges[3])] == NULL)
1229  knotVectors[KnotInd(edges[3])] =
1230  new KnotVector(*patches[p]->GetKV(1));
1231  if (knotVectors[KnotInd(edges[8])] == NULL)
1232  knotVectors[KnotInd(edges[8])] =
1233  new KnotVector(*patches[p]->GetKV(2));
1234  }
1235  }
1236  Order = knotVectors[0]->GetOrder();
1237  }
1238
1239  GenerateOffsets();
1240  CountElements();
1241  CountBdrElements();
1242  // NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs
1243
1244  skip_comment_lines(input, '#');
1245
1246  // Check for a list of mesh elements
1247  if (patches.Size() == 0)
1248  {
1249  input >> ws >> ident;
1250  }
1251  if (patches.Size() == 0 && ident == "mesh_elements")
1252  {
1253  input >> NumOfActiveElems;
1254  activeElem.SetSize(GetGNE());
1255  activeElem = false;
1256  int glob_elem;
1257  for (int i = 0; i < NumOfActiveElems; i++)
1258  {
1259  input >> glob_elem;
1260  activeElem[glob_elem] = true;
1261  }
1262
1263  skip_comment_lines(input, '#');
1264  input >> ws >> ident;
1265  }
1266  else
1267  {
1268  NumOfActiveElems = NumOfElements;
1269  activeElem.SetSize(NumOfElements);
1270  activeElem = true;
1271  }
1272
1273  GenerateActiveVertices();
1274  GenerateElementDofTable();
1275  GenerateActiveBdrElems();
1276  GenerateBdrElementDofTable();
1277
1278  if (patches.Size() == 0)
1279  {
1280  // weights
1281  if (ident == "weights")
1282  {
1284  }
1285  else // e.g. ident = "unitweights" or "autoweights"
1286  {
1287  weights.SetSize(GetNDof());
1288  weights = 1.0;
1289  }
1290  }
1291 }
1292
1294 {
1295  Order = Order_;
1296
1297  patchTopo = parent->patchTopo;
1298  own_topo = 0;
1299
1300  parent->edge_to_knot.Copy(edge_to_knot);
1301
1302  NumOfKnotVectors = parent->GetNKV();
1303  knotVectors.SetSize(NumOfKnotVectors);
1304  for (int i = 0; i < NumOfKnotVectors; i++)
1305  {
1306  knotVectors[i] =
1307  parent->GetKnotVector(i)->DegreeElevate(Order-parent->GetOrder());
1308  }
1309
1310  // copy some data from parent
1311  NumOfElements = parent->NumOfElements;
1312  NumOfBdrElements = parent->NumOfBdrElements;
1313
1314  GenerateOffsets(); // dof offsets will be different from parent
1315
1316  NumOfActiveVertices = parent->NumOfActiveVertices;
1317  NumOfActiveElems = parent->NumOfActiveElems;
1318  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1319  parent->activeVert.Copy(activeVert);
1320  parent->activeElem.Copy(activeElem);
1321  parent->activeBdrElem.Copy(activeBdrElem);
1322
1323  GenerateElementDofTable();
1324  GenerateBdrElementDofTable();
1325
1326  weights.SetSize(GetNDof());
1327  weights = 1.0;
1328 }
1329
1330 NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
1331 {
1332  NURBSExtension *parent = mesh_array[0]->NURBSext;
1333
1334  if (!parent->own_topo)
1335  mfem_error("NURBSExtension::NURBSExtension :\n"
1336  " parent does not own the patch topology!");
1337  patchTopo = parent->patchTopo;
1338  own_topo = 1;
1339  parent->own_topo = 0;
1340
1341  parent->edge_to_knot.Copy(edge_to_knot);
1342
1343  Order = parent->GetOrder();
1344
1345  NumOfKnotVectors = parent->GetNKV();
1346  knotVectors.SetSize(NumOfKnotVectors);
1347  for (int i = 0; i < NumOfKnotVectors; i++)
1348  {
1349  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1350  }
1351
1352  GenerateOffsets();
1353  CountElements();
1354  CountBdrElements();
1355
1356  // assuming the meshes define a partitioning of all the elements
1357  NumOfActiveElems = NumOfElements;
1358  activeElem.SetSize(NumOfElements);
1359  activeElem = true;
1360
1361  GenerateActiveVertices();
1362  GenerateElementDofTable();
1363  GenerateActiveBdrElems();
1364  GenerateBdrElementDofTable();
1365
1366  weights.SetSize(GetNDof());
1367  MergeWeights(mesh_array, num_pieces);
1368 }
1369
1371 {
1372  if (patches.Size() == 0)
1373  {
1374  delete bel_dof;
1375  delete el_dof;
1376  }
1377
1378  for (int i = 0; i < knotVectors.Size(); i++)
1379  {
1380  delete knotVectors[i];
1381  }
1382
1383  for (int i = 0; i < patches.Size(); i++)
1384  {
1385  delete patches[i];
1386  }
1387
1388  if (own_topo)
1389  {
1390  delete patchTopo;
1391  }
1392 }
1393
1394 void NURBSExtension::Print(std::ostream &out) const
1395 {
1396  patchTopo->PrintTopo(out, edge_to_knot);
1397  if (patches.Size() == 0)
1398  {
1399  out << "\nknotvectors\n" << NumOfKnotVectors << '\n';
1400  for (int i = 0; i < NumOfKnotVectors; i++)
1401  {
1402  knotVectors[i]->Print(out);
1403  }
1404
1405  if (NumOfActiveElems < NumOfElements)
1406  {
1407  out << "\nmesh_elements\n" << NumOfActiveElems << '\n';
1408  for (int i = 0; i < NumOfElements; i++)
1409  if (activeElem[i])
1410  {
1411  out << i << '\n';
1412  }
1413  }
1414
1415  out << "\nweights\n";
1416  weights.Print(out, 1);
1417  }
1418  else
1419  {
1420  out << "\npatches\n";
1421  for (int p = 0; p < patches.Size(); p++)
1422  {
1423  out << "\n# patch " << p << "\n\n";
1424  patches[p]->Print(out);
1425  }
1426  }
1427 }
1428
1430 {
1431  out <<
1432  "NURBS Mesh entity sizes:\n"
1433  "Dimension = " << Dimension() << "\n"
1434  "Order = " << GetOrder() << "\n"
1435  "NumOfKnotVectors = " << GetNKV() << "\n"
1436  "NumOfPatches = " << GetNP() << "\n"
1437  "NumOfBdrPatches = " << GetNBP() << "\n"
1438  "NumOfVertices = " << GetGNV() << "\n"
1439  "NumOfElements = " << GetGNE() << "\n"
1440  "NumOfBdrElements = " << GetGNBE() << "\n"
1441  "NumOfDofs = " << GetNTotalDof() << "\n"
1442  "NumOfActiveVertices = " << GetNV() << "\n"
1443  "NumOfActiveElems = " << GetNE() << "\n"
1444  "NumOfActiveBdrElems = " << GetNBE() << "\n"
1445  "NumOfActiveDofs = " << GetNDof() << '\n';
1446  for (int i = 0; i < NumOfKnotVectors; i++)
1447  {
1448  out << ' ' << i + 1 << ") ";
1449  knotVectors[i]->Print(out);
1450  }
1451  cout << endl;
1452 }
1453
1455 {
1456  int vert[8], nv, g_el, nx, ny, nz, dim = Dimension();
1457
1458  NURBSPatchMap p2g(this);
1459  KnotVector *kv[3];
1460
1461  g_el = 0;
1462  activeVert.SetSize(GetGNV());
1463  activeVert = -1;
1464  for (int p = 0; p < GetNP(); p++)
1465  {
1466  p2g.SetPatchVertexMap(p, kv);
1467
1468  nx = p2g.nx();
1469  ny = p2g.ny();
1470  nz = (dim == 3) ? p2g.nz() : 1;
1471
1472  for (int k = 0; k < nz; k++)
1473  {
1474  for (int j = 0; j < ny; j++)
1475  {
1476  for (int i = 0; i < nx; i++)
1477  {
1478  if (activeElem[g_el])
1479  {
1480  if (dim == 2)
1481  {
1482  vert[0] = p2g(i, j );
1483  vert[1] = p2g(i+1,j );
1484  vert[2] = p2g(i+1,j+1);
1485  vert[3] = p2g(i, j+1);
1486  nv = 4;
1487  }
1488  else
1489  {
1490  vert[0] = p2g(i, j, k);
1491  vert[1] = p2g(i+1,j, k);
1492  vert[2] = p2g(i+1,j+1,k);
1493  vert[3] = p2g(i, j+1,k);
1494
1495  vert[4] = p2g(i, j, k+1);
1496  vert[5] = p2g(i+1,j, k+1);
1497  vert[6] = p2g(i+1,j+1,k+1);
1498  vert[7] = p2g(i, j+1,k+1);
1499  nv = 8;
1500  }
1501
1502  for (int v = 0; v < nv; v++)
1503  {
1504  activeVert[vert[v]] = 1;
1505  }
1506  }
1507  g_el++;
1508  }
1509  }
1510  }
1511  }
1512
1513  NumOfActiveVertices = 0;
1514  for (int i = 0; i < GetGNV(); i++)
1515  if (activeVert[i] == 1)
1516  {
1517  activeVert[i] = NumOfActiveVertices++;
1518  }
1519 }
1520
1522 {
1523  int dim = Dimension();
1524  Array<KnotVector *> kv(dim);
1525
1526  activeBdrElem.SetSize(GetGNBE());
1527  if (GetGNE() == GetNE())
1528  {
1529  activeBdrElem = true;
1530  NumOfActiveBdrElems = GetGNBE();
1531  return;
1532  }
1533  activeBdrElem = false;
1534  NumOfActiveBdrElems = 0;
1535  // the mesh will generate the actual boundary including boundary
1536  // elements that are not on boundary patches. we use this for
1537  // visualization of processor boundaries
1538
1539  // TODO: generate actual boundary?
1540 }
1541
1542
1543 void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces)
1544 {
1545  Array<int> lelem_elem;
1546
1547  for (int i = 0; i < num_pieces; i++)
1548  {
1549  NURBSExtension *lext = mesh_array[i]->NURBSext;
1550
1551  lext->GetElementLocalToGlobal(lelem_elem);
1552
1553  for (int lel = 0; lel < lext->GetNE(); lel++)
1554  {
1555  int gel = lelem_elem[lel];
1556
1557  int nd = el_dof->RowSize(gel);
1558  int *gdofs = el_dof->GetRow(gel);
1559  int *ldofs = lext->el_dof->GetRow(lel);
1560  for (int j = 0; j < nd; j++)
1561  {
1562  weights(gdofs[j]) = lext->weights(ldofs[j]);
1563  }
1564  }
1565  }
1566 }
1567
1569  GridFunction *gf_array[], int num_pieces, GridFunction &merged)
1570 {
1571  FiniteElementSpace *gfes = merged.FESpace();
1572  Array<int> lelem_elem, dofs;
1573  Vector lvec;
1574
1575  for (int i = 0; i < num_pieces; i++)
1576  {
1577  FiniteElementSpace *lfes = gf_array[i]->FESpace();
1578  NURBSExtension *lext = lfes->GetMesh()->NURBSext;
1579
1580  lext->GetElementLocalToGlobal(lelem_elem);
1581
1582  for (int lel = 0; lel < lext->GetNE(); lel++)
1583  {
1584  lfes->GetElementVDofs(lel, dofs);
1585  gf_array[i]->GetSubVector(dofs, lvec);
1586
1587  gfes->GetElementVDofs(lelem_elem[lel], dofs);
1588  merged.SetSubVector(dofs, lvec);
1589  }
1590  }
1591 }
1592
1594 {
1595  Array<int> edges;
1596  Array<int> oedge;
1597
1598  for (int p = 0; p < GetNP(); p++)
1599  {
1600  patchTopo->GetElementEdges(p, edges, oedge);
1601
1602  for (int i = 0; i < edges.Size(); i++)
1603  {
1604  edges[i] = edge_to_knot[edges[i]];
1605  if (oedge[i] < 0)
1606  {
1607  edges[i] = -1 - edges[i];
1608  }
1609  }
1610
1611  if ((Dimension() == 2 &&
1612  (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
1613
1614  (Dimension() == 3 &&
1615  (edges[0] != edges[2] || edges[0] != edges[4] ||
1616  edges[0] != edges[6] || edges[1] != edges[3] ||
1617  edges[1] != edges[5] || edges[1] != edges[7] ||
1618  edges[8] != edges[9] || edges[8] != edges[10] ||
1619  edges[8] != edges[11])))
1620  {
1621  cerr << "NURBSExtension::CheckPatch (patch = " << p
1622  << ")\n Inconsistent edge-to-knot mapping!\n";
1623  mfem_error();
1624  }
1625
1626  if ((Dimension() == 2 &&
1627  (edges[0] < 0 || edges[1] < 0)) ||
1628
1629  (Dimension() == 3 &&
1630  (edges[0] < 0 || edges[3] < 0 || edges[8] < 0)))
1631  {
1632  cerr << "NURBSExtension::CheckPatch (patch = " << p
1633  << ") : Bad orientation!\n";
1634  mfem_error();
1635  }
1636  }
1637 }
1638
1640 {
1641  Array<int> edges;
1642  Array<int> oedge;
1643
1644  for (int p = 0; p < GetNBP(); p++)
1645  {
1646  patchTopo->GetBdrElementEdges(p, edges, oedge);
1647
1648  for (int i = 0; i < edges.Size(); i++)
1649  {
1650  edges[i] = edge_to_knot[edges[i]];
1651  if (oedge[i] < 0)
1652  {
1653  edges[i] = -1 - edges[i];
1654  }
1655  }
1656
1657  if ((Dimension() == 2 && (edges[0] < 0)) ||
1658  (Dimension() == 3 && (edges[0] < 0 || edges[1] < 0)))
1659  {
1660  cerr << "NURBSExtension::CheckBdrPatch (boundary patch = "
1661  << p << ") : Bad orientation!\n";
1662  mfem_error();
1663  }
1664  }
1665 }
1666
1668 {
1669  Array<int> edges;
1670  Array<int> orient;
1671
1672  kv.SetSize(Dimension());
1673  patchTopo->GetElementEdges(p, edges, orient);
1674
1675  if (Dimension() == 2)
1676  {
1677  kv[0] = KnotVec(edges[0]);
1678  kv[1] = KnotVec(edges[1]);
1679  }
1680  else
1681  {
1682  kv[0] = KnotVec(edges[0]);
1683  kv[1] = KnotVec(edges[3]);
1684  kv[2] = KnotVec(edges[8]);
1685  }
1686 }
1687
1689 {
1690  Array<int> edges;
1691  Array<int> orient;
1692
1693  kv.SetSize(Dimension() - 1);
1694  patchTopo->GetBdrElementEdges(p, edges, orient);
1695
1696  if (Dimension() == 2)
1697  {
1698  kv[0] = KnotVec(edges[0]);
1699  }
1700  else
1701  {
1702  kv[0] = KnotVec(edges[0]);
1703  kv[1] = KnotVec(edges[1]);
1704  }
1705 }
1706
1708 {
1709  int nv = patchTopo->GetNV();
1710  int ne = patchTopo->GetNEdges();
1711  int nf = patchTopo->GetNFaces();
1712  int np = patchTopo->GetNE();
1713  int meshCounter, spaceCounter, dim = Dimension();
1714
1715  Array<int> edges;
1716  Array<int> orient;
1717
1718  v_meshOffsets.SetSize(nv);
1719  e_meshOffsets.SetSize(ne);
1720  f_meshOffsets.SetSize(nf);
1721  p_meshOffsets.SetSize(np);
1722
1723  v_spaceOffsets.SetSize(nv);
1724  e_spaceOffsets.SetSize(ne);
1725  f_spaceOffsets.SetSize(nf);
1726  p_spaceOffsets.SetSize(np);
1727
1728  // Get vertex offsets
1729  for (meshCounter = 0; meshCounter < nv; meshCounter++)
1730  {
1731  v_meshOffsets[meshCounter] = meshCounter;
1732  v_spaceOffsets[meshCounter] = meshCounter;
1733  }
1734  spaceCounter = meshCounter;
1735
1736  // Get edge offsets
1737  for (int e = 0; e < ne; e++)
1738  {
1739  e_meshOffsets[e] = meshCounter;
1740  e_spaceOffsets[e] = spaceCounter;
1741  meshCounter += KnotVec(e)->GetNE() - 1;
1742  spaceCounter += KnotVec(e)->GetNCP() - 2;
1743  }
1744
1745  // Get face offsets
1746  for (int f = 0; f < nf; f++)
1747  {
1748  f_meshOffsets[f] = meshCounter;
1749  f_spaceOffsets[f] = spaceCounter;
1750
1751  patchTopo->GetFaceEdges(f, edges, orient);
1752
1753  meshCounter +=
1754  (KnotVec(edges[0])->GetNE() - 1) *
1755  (KnotVec(edges[1])->GetNE() - 1);
1756  spaceCounter +=
1757  (KnotVec(edges[0])->GetNCP() - 2) *
1758  (KnotVec(edges[1])->GetNCP() - 2);
1759  }
1760
1761  // Get patch offsets
1762  for (int p = 0; p < np; p++)
1763  {
1764  p_meshOffsets[p] = meshCounter;
1765  p_spaceOffsets[p] = spaceCounter;
1766
1767  patchTopo->GetElementEdges(p, edges, orient);
1768
1769  if (dim == 2)
1770  {
1771  meshCounter +=
1772  (KnotVec(edges[0])->GetNE() - 1) *
1773  (KnotVec(edges[1])->GetNE() - 1);
1774  spaceCounter +=
1775  (KnotVec(edges[0])->GetNCP() - 2) *
1776  (KnotVec(edges[1])->GetNCP() - 2);
1777  }
1778  else
1779  {
1780  meshCounter +=
1781  (KnotVec(edges[0])->GetNE() - 1) *
1782  (KnotVec(edges[3])->GetNE() - 1) *
1783  (KnotVec(edges[8])->GetNE() - 1);
1784  spaceCounter +=
1785  (KnotVec(edges[0])->GetNCP() - 2) *
1786  (KnotVec(edges[3])->GetNCP() - 2) *
1787  (KnotVec(edges[8])->GetNCP() - 2);
1788  }
1789  }
1790  NumOfVertices = meshCounter;
1791  NumOfDofs = spaceCounter;
1792 }
1793
1795 {
1796  int dim = Dimension();
1797  Array<KnotVector *> kv(dim);
1798
1799  NumOfElements = 0;
1800  for (int p = 0; p < GetNP(); p++)
1801  {
1802  GetPatchKnotVectors(p, kv);
1803
1804  int ne = kv[0]->GetNE();
1805  for (int d = 1; d < dim; d++)
1806  {
1807  ne *= kv[d]->GetNE();
1808  }
1809
1810  NumOfElements += ne;
1811  }
1812 }
1813
1815 {
1816  int dim = Dimension() - 1;
1817  Array<KnotVector *> kv(dim);
1818
1819  NumOfBdrElements = 0;
1820  for (int p = 0; p < GetNBP(); p++)
1821  {
1822  GetBdrPatchKnotVectors(p, kv);
1823
1824  int ne = kv[0]->GetNE();
1825  for (int d = 1; d < dim; d++)
1826  {
1827  ne *= kv[d]->GetNE();
1828  }
1829
1830  NumOfBdrElements += ne;
1831  }
1832 }
1833
1835 {
1836  elements.SetSize(GetNE());
1837
1838  if (Dimension() == 2)
1839  {
1840  Get2DElementTopo(elements);
1841  }
1842  else
1843  {
1844  Get3DElementTopo(elements);
1845  }
1846 }
1847
1849 {
1850  int el = 0;
1851  int eg = 0;
1852  int ind[4];
1853  NURBSPatchMap p2g(this);
1854  KnotVector *kv[2];
1855
1856  for (int p = 0; p < GetNP(); p++)
1857  {
1858  p2g.SetPatchVertexMap(p, kv);
1859  int nx = p2g.nx();
1860  int ny = p2g.ny();
1861
1862  int patch_attr = patchTopo->GetAttribute(p);
1863
1864  for (int j = 0; j < ny; j++)
1865  {
1866  for (int i = 0; i < nx; i++)
1867  {
1868  if (activeElem[eg])
1869  {
1870  ind[0] = activeVert[p2g(i, j )];
1871  ind[1] = activeVert[p2g(i+1,j )];
1872  ind[2] = activeVert[p2g(i+1,j+1)];
1873  ind[3] = activeVert[p2g(i, j+1)];
1874
1875  elements[el] = new Quadrilateral(ind, patch_attr);
1876  el++;
1877  }
1878  eg++;
1879  }
1880  }
1881  }
1882 }
1883
1885 {
1886  int el = 0;
1887  int eg = 0;
1888  int ind[8];
1889  NURBSPatchMap p2g(this);
1890  KnotVector *kv[3];
1891
1892  for (int p = 0; p < GetNP(); p++)
1893  {
1894  p2g.SetPatchVertexMap(p, kv);
1895  int nx = p2g.nx();
1896  int ny = p2g.ny();
1897  int nz = p2g.nz();
1898
1899  int patch_attr = patchTopo->GetAttribute(p);
1900
1901  for (int k = 0; k < nz; k++)
1902  {
1903  for (int j = 0; j < ny; j++)
1904  {
1905  for (int i = 0; i < nx; i++)
1906  {
1907  if (activeElem[eg])
1908  {
1909  ind[0] = activeVert[p2g(i, j, k)];
1910  ind[1] = activeVert[p2g(i+1,j, k)];
1911  ind[2] = activeVert[p2g(i+1,j+1,k)];
1912  ind[3] = activeVert[p2g(i, j+1,k)];
1913
1914  ind[4] = activeVert[p2g(i, j, k+1)];
1915  ind[5] = activeVert[p2g(i+1,j, k+1)];
1916  ind[6] = activeVert[p2g(i+1,j+1,k+1)];
1917  ind[7] = activeVert[p2g(i, j+1,k+1)];
1918
1919  elements[el] = new Hexahedron(ind, patch_attr);
1920  el++;
1921  }
1922  eg++;
1923  }
1924  }
1925  }
1926  }
1927 }
1928
1930 {
1931  boundary.SetSize(GetNBE());
1932
1933  if (Dimension() == 2)
1934  {
1935  Get2DBdrElementTopo(boundary);
1936  }
1937  else
1938  {
1939  Get3DBdrElementTopo(boundary);
1940  }
1941 }
1942
1944 {
1945  int g_be, l_be;
1946  int ind[2], okv[1];
1947  NURBSPatchMap p2g(this);
1948  KnotVector *kv[1];
1949
1950  g_be = l_be = 0;
1951  for (int b = 0; b < GetNBP(); b++)
1952  {
1953  p2g.SetBdrPatchVertexMap(b, kv, okv);
1954  int nx = p2g.nx();
1955
1956  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
1957
1958  for (int i = 0; i < nx; i++)
1959  {
1960  if (activeBdrElem[g_be])
1961  {
1962  int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
1963  ind[0] = activeVert[p2g[_i ]];
1964  ind[1] = activeVert[p2g[_i+1]];
1965
1966  boundary[l_be] = new Segment(ind, bdr_patch_attr);
1967  l_be++;
1968  }
1969  g_be++;
1970  }
1971  }
1972 }
1973
1975 {
1976  int g_be, l_be;
1977  int ind[4], okv[2];
1978  NURBSPatchMap p2g(this);
1979  KnotVector *kv[2];
1980
1981  g_be = l_be = 0;
1982  for (int b = 0; b < GetNBP(); b++)
1983  {
1984  p2g.SetBdrPatchVertexMap(b, kv, okv);
1985  int nx = p2g.nx();
1986  int ny = p2g.ny();
1987
1988  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
1989
1990  for (int j = 0; j < ny; j++)
1991  {
1992  int _j = (okv[1] >= 0) ? j : (ny - 1 - j);
1993  for (int i = 0; i < nx; i++)
1994  {
1995  if (activeBdrElem[g_be])
1996  {
1997  int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
1998  ind[0] = activeVert[p2g(_i, _j )];
1999  ind[1] = activeVert[p2g(_i+1,_j )];
2000  ind[2] = activeVert[p2g(_i+1,_j+1)];
2001  ind[3] = activeVert[p2g(_i, _j+1)];
2002
2003  boundary[l_be] = new Quadrilateral(ind, bdr_patch_attr);
2004  l_be++;
2005  }
2006  g_be++;
2007  }
2008  }
2009  }
2010 }
2011
2013 {
2014  activeDof.SetSize(GetNTotalDof());
2015  activeDof = 0;
2016
2017  if (Dimension() == 2)
2018  {
2019  Generate2DElementDofTable();
2020  }
2021  else
2022  {
2023  Generate3DElementDofTable();
2024  }
2025
2026  NumOfActiveDofs = 0;
2027  for (int d = 0; d < GetNTotalDof(); d++)
2028  if (activeDof[d])
2029  {
2030  NumOfActiveDofs++;
2031  activeDof[d] = NumOfActiveDofs;
2032  }
2033
2034  int *dof = el_dof->GetJ();
2035  int ndof = el_dof->Size_of_connections();
2036  for (int i = 0; i < ndof; i++)
2037  {
2038  dof[i] = activeDof[dof[i]] - 1;
2039  }
2040 }
2041
2043 {
2044  int el = 0;
2045  int eg = 0;
2046  KnotVector *kv[2];
2047  NURBSPatchMap p2g(this);
2048
2049  el_dof = new Table(NumOfActiveElems, (Order + 1)*(Order + 1));
2050  el_to_patch.SetSize(NumOfActiveElems);
2051  el_to_IJK.SetSize(NumOfActiveElems, 2);
2052
2053  for (int p = 0; p < GetNP(); p++)
2054  {
2055  p2g.SetPatchDofMap(p, kv);
2056
2058  for (int j = 0; j < kv[1]->GetNKS(); j++)
2059  {
2060  if (kv[1]->isElement(j))
2061  {
2062  for (int i = 0; i < kv[0]->GetNKS(); i++)
2063  {
2064  if (kv[0]->isElement(i))
2065  {
2066  if (activeElem[eg])
2067  {
2068  int *dofs = el_dof->GetRow(el);
2069  int idx = 0;
2070  for (int jj = 0; jj <= Order; jj++)
2071  {
2072  for (int ii = 0; ii <= Order; ii++)
2073  {
2074  dofs[idx] = p2g(i+ii,j+jj);
2075  activeDof[dofs[idx]] = 1;
2076  idx++;
2077  }
2078  }
2079  el_to_patch[el] = p;
2080  el_to_IJK(el,0) = i;
2081  el_to_IJK(el,1) = j;
2082
2083  el++;
2084  }
2085  eg++;
2086  }
2087  }
2088  }
2089  }
2090  }
2091 }
2092
2094 {
2095  int el = 0;
2096  int eg = 0;
2097  int idx;
2098  KnotVector *kv[3];
2099  NURBSPatchMap p2g(this);
2100
2101  el_dof = new Table(NumOfActiveElems, (Order + 1)*(Order + 1)*(Order + 1));
2102  el_to_patch.SetSize(NumOfActiveElems);
2103  el_to_IJK.SetSize(NumOfActiveElems, 3);
2104
2105  for (int p = 0; p < GetNP(); p++)
2106  {
2107  p2g.SetPatchDofMap(p, kv);
2108
2110  for (int k = 0; k < kv[2]->GetNKS(); k++)
2111  {
2112  if (kv[2]->isElement(k))
2113  {
2114  for (int j = 0; j < kv[1]->GetNKS(); j++)
2115  {
2116  if (kv[1]->isElement(j))
2117  {
2118  for (int i = 0; i < kv[0]->GetNKS(); i++)
2119  {
2120  if (kv[0]->isElement(i))
2121  {
2122  if (activeElem[eg])
2123  {
2124  idx = 0;
2125  int *dofs = el_dof->GetRow(el);
2126  for (int kk = 0; kk <= Order; kk++)
2127  {
2128  for (int jj = 0; jj <= Order; jj++)
2129  {
2130  for (int ii = 0; ii <= Order; ii++)
2131  {
2132  dofs[idx] = p2g(i+ii,j+jj,k+kk);
2133  activeDof[dofs[idx]] = 1;
2134  idx++;
2135  }
2136  }
2137  }
2138
2139  el_to_patch[el] = p;
2140  el_to_IJK(el,0) = i;
2141  el_to_IJK(el,1) = j;
2142  el_to_IJK(el,2) = k;
2143
2144  el++;
2145  }
2146  eg++;
2147  }
2148  }
2149  }
2150  }
2151  }
2152  }
2153  }
2154 }
2155
2157 {
2158  if (Dimension() == 2)
2159  {
2160  Generate2DBdrElementDofTable();
2161  }
2162  else
2163  {
2164  Generate3DBdrElementDofTable();
2165  }
2166
2167  int *dof = bel_dof->GetJ();
2168  int ndof = bel_dof->Size_of_connections();
2169  for (int i = 0; i < ndof; i++)
2170  {
2171  dof[i] = activeDof[dof[i]] - 1;
2172  }
2173 }
2174
2176 {
2177  int gbe = 0;
2178  int lbe = 0, okv[1];
2179  KnotVector *kv[1];
2180  NURBSPatchMap p2g(this);
2181
2182  bel_dof = new Table(NumOfActiveBdrElems, Order + 1);
2183  bel_to_patch.SetSize(NumOfActiveBdrElems);
2184  bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
2185
2186  for (int b = 0; b < GetNBP(); b++)
2187  {
2188  p2g.SetBdrPatchDofMap(b, kv, okv);
2189  int nx = p2g.nx(); // NCP-1
2190
2192  int nks0 = kv[0]->GetNKS();
2193  for (int i = 0; i < nks0; i++)
2194  {
2195  if (kv[0]->isElement(i))
2196  {
2197  if (activeBdrElem[gbe])
2198  {
2199  int *dofs = bel_dof->GetRow(lbe);
2200  int idx = 0;
2201  for (int ii = 0; ii <= Order; ii++)
2202  {
2203  dofs[idx] = p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)];
2204  idx++;
2205  }
2206  bel_to_patch[lbe] = b;
2207  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2208  lbe++;
2209  }
2210  gbe++;
2211  }
2212  }
2213  }
2214 }
2215
2217 {
2218  int gbe = 0;
2219  int lbe = 0, okv[2];
2220  KnotVector *kv[2];
2221  NURBSPatchMap p2g(this);
2222
2223  bel_dof = new Table(NumOfActiveBdrElems, (Order + 1)*(Order + 1));
2224  bel_to_patch.SetSize(NumOfActiveBdrElems);
2225  bel_to_IJK.SetSize(NumOfActiveBdrElems, 2);
2226
2227  for (int b = 0; b < GetNBP(); b++)
2228  {
2229  p2g.SetBdrPatchDofMap(b, kv, okv);
2230  int nx = p2g.nx(); // NCP0-1
2231  int ny = p2g.ny(); // NCP1-1
2232
2234  int nks0 = kv[0]->GetNKS();
2235  int nks1 = kv[1]->GetNKS();
2236  for (int j = 0; j < nks1; j++)
2237  {
2238  if (kv[1]->isElement(j))
2239  {
2240  for (int i = 0; i < nks0; i++)
2241  {
2242  if (kv[0]->isElement(i))
2243  {
2244  if (activeBdrElem[gbe])
2245  {
2246  int *dofs = bel_dof->GetRow(lbe);
2247  int idx = 0;
2248  for (int jj = 0; jj <= Order; jj++)
2249  {
2250  int _jj = (okv[1] >= 0) ? (j+jj) : (ny-j-jj);
2251  for (int ii = 0; ii <= Order; ii++)
2252  {
2253  int _ii = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
2254  dofs[idx] = p2g(_ii,_jj);
2255  idx++;
2256  }
2257  }
2258  bel_to_patch[lbe] = b;
2259  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2260  bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
2261  lbe++;
2262  }
2263  gbe++;
2264  }
2265  }
2266  }
2267  }
2268  }
2269 }
2270
2272 {
2273  lvert_vert.SetSize(GetNV());
2274  for (int gv = 0; gv < GetGNV(); gv++)
2275  if (activeVert[gv] >= 0)
2276  {
2277  lvert_vert[activeVert[gv]] = gv;
2278  }
2279 }
2280
2282 {
2283  lelem_elem.SetSize(GetNE());
2284  for (int le = 0, ge = 0; ge < GetGNE(); ge++)
2285  if (activeElem[ge])
2286  {
2287  lelem_elem[le++] = ge;
2288  }
2289 }
2290
2292 {
2293  const NURBSFiniteElement *NURBSFE =
2294  dynamic_cast<const NURBSFiniteElement *>(FE);
2295
2296  if (NURBSFE->GetElement() != i)
2297  {
2298  Array<int> dofs;
2299  NURBSFE->SetIJK(el_to_IJK.GetRow(i));
2300  if (el_to_patch[i] != NURBSFE->GetPatch())
2301  {
2302  GetPatchKnotVectors(el_to_patch[i], NURBSFE->KnotVectors());
2303  NURBSFE->SetPatch(el_to_patch[i]);
2304  }
2305  el_dof->GetRow(i, dofs);
2306  weights.GetSubVector(dofs, NURBSFE->Weights());
2307  NURBSFE->SetElement(i);
2308  }
2309 }
2310
2312 {
2313  const NURBSFiniteElement *NURBSFE =
2314  dynamic_cast<const NURBSFiniteElement *>(BE);
2315
2316  if (NURBSFE->GetElement() != i)
2317  {
2318  Array<int> dofs;
2319  NURBSFE->SetIJK(bel_to_IJK.GetRow(i));
2320  if (bel_to_patch[i] != NURBSFE->GetPatch())
2321  {
2322  GetBdrPatchKnotVectors(bel_to_patch[i], NURBSFE->KnotVectors());
2323  NURBSFE->SetPatch(bel_to_patch[i]);
2324  }
2325  bel_dof->GetRow(i, dofs);
2326  weights.GetSubVector(dofs, NURBSFE->Weights());
2327  NURBSFE->SetElement(i);
2328  }
2329 }
2330
2332 {
2333  delete el_dof;
2334  delete bel_dof;
2335
2336  if (patches.Size() == 0)
2337  {
2338  GetPatchNets(Nodes);
2339  }
2340 }
2341
2343 {
2344  if (patches.Size() == 0) { return; }
2345
2346  SetSolutionVector(Nodes);
2347  patches.SetSize(0);
2348 }
2349
2351 {
2352  if (patches.Size() == 0)
2353  mfem_error("NURBSExtension::SetKnotsFromPatches :"
2354  " No patches available!");
2355
2357
2358  for (int p = 0; p < patches.Size(); p++)
2359  {
2360  GetPatchKnotVectors(p, kv);
2361
2362  for (int i = 0; i < kv.Size(); i++)
2363  {
2364  *kv[i] = *patches[p]->GetKV(i);
2365  }
2366  }
2367
2368  Order = knotVectors[0]->GetOrder();
2369  for (int i = 0; i < NumOfKnotVectors; i++)
2370  {
2371  if (Order != knotVectors[i]->GetOrder())
2372  mfem_error("NURBSExtension::Reset :\n"
2373  " Variable orders are not supported!");
2374  }
2375
2376  GenerateOffsets();
2377  CountElements();
2378  CountBdrElements();
2379
2380  // all elements must be active
2381  NumOfActiveElems = NumOfElements;
2382  activeElem.SetSize(NumOfElements);
2383  activeElem = true;
2384
2385  GenerateActiveVertices();
2386  GenerateElementDofTable();
2387  GenerateActiveBdrElems();
2388  GenerateBdrElementDofTable();
2389 }
2390
2392 {
2393  for (int p = 0; p < patches.Size(); p++)
2394  {
2395  patches[p]->DegreeElevate(t);
2396  }
2397 }
2398
2400 {
2401  for (int p = 0; p < patches.Size(); p++)
2402  {
2403  patches[p]->UniformRefinement();
2404  }
2405 }
2406
2408 {
2409  Array<int> edges;
2410  Array<int> orient;
2411
2412  Array<KnotVector *> pkv(Dimension());
2413
2414  for (int p = 0; p < patches.Size(); p++)
2415  {
2416  patchTopo->GetElementEdges(p, edges, orient);
2417
2418  if (Dimension()==2)
2419  {
2420  pkv[0] = kv[KnotInd(edges[0])];
2421  pkv[1] = kv[KnotInd(edges[1])];
2422  }
2423  else
2424  {
2425  pkv[0] = kv[KnotInd(edges[0])];
2426  pkv[1] = kv[KnotInd(edges[3])];
2427  pkv[2] = kv[KnotInd(edges[8])];
2428  }
2429
2430  patches[p]->KnotInsert(pkv);
2431  }
2432 }
2433
2435 {
2436  if (Dimension() == 2)
2437  {
2438  Get2DPatchNets(coords);
2439  }
2440  else
2441  {
2442  Get3DPatchNets(coords);
2443  }
2444 }
2445
2447 {
2448  Array<KnotVector *> kv(2);
2449  NURBSPatchMap p2g(this);
2450
2451  patches.SetSize(GetNP());
2452  for (int p = 0; p < GetNP(); p++)
2453  {
2454  p2g.SetPatchDofMap(p, kv);
2455  patches[p] = new NURBSPatch(kv, 2+1);
2456  NURBSPatch &Patch = *patches[p];
2457
2458  for (int j = 0; j < kv[1]->GetNCP(); j++)
2459  {
2460  for (int i = 0; i < kv[0]->GetNCP(); i++)
2461  {
2462  int l = p2g(i,j);
2463  for (int d = 0; d < 2; d++)
2464  {
2465  Patch(i,j,d) = coords(l*2 + d)*weights(l);
2466  }
2467  Patch(i,j,2) = weights(l);
2468  }
2469  }
2470  }
2471 }
2472
2474 {
2475  Array<KnotVector *> kv(3);
2476  NURBSPatchMap p2g(this);
2477
2478  patches.SetSize(GetNP());
2479  for (int p = 0; p < GetNP(); p++)
2480  {
2481  p2g.SetPatchDofMap(p, kv);
2482  patches[p] = new NURBSPatch(kv, 3+1);
2483  NURBSPatch &Patch = *patches[p];
2484
2485  for (int k = 0; k < kv[2]->GetNCP(); k++)
2486  {
2487  for (int j = 0; j < kv[1]->GetNCP(); j++)
2488  {
2489  for (int i = 0; i < kv[0]->GetNCP(); i++)
2490  {
2491  int l = p2g(i,j,k);
2492  for (int d = 0; d < 3; d++)
2493  {
2494  Patch(i,j,k,d) = coords(l*3 + d)*weights(l);
2495  }
2496  Patch(i,j,k,3) = weights(l);
2497  }
2498  }
2499  }
2500  }
2501 }
2502
2504 {
2505  if (Dimension() == 2)
2506  {
2507  Set2DSolutionVector(coords);
2508  }
2509  else
2510  {
2511  Set3DSolutionVector(coords);
2512  }
2513 }
2514
2516 {
2517  Array<KnotVector *> kv(2);
2518  NURBSPatchMap p2g(this);
2519
2520  weights.SetSize(GetNDof());
2521  for (int p = 0; p < GetNP(); p++)
2522  {
2523  p2g.SetPatchDofMap(p, kv);
2524  NURBSPatch &Patch = *patches[p];
2525
2526  for (int j = 0; j < kv[1]->GetNCP(); j++)
2527  {
2528  for (int i = 0; i < kv[0]->GetNCP(); i++)
2529  {
2530  int l = p2g(i,j);
2531  for (int d = 0; d < 2; d++)
2532  {
2533  coords(l*2 + d) = Patch(i,j,d)/Patch(i,j,2);
2534  }
2535  weights(l) = Patch(i,j,2);
2536  }
2537  }
2538  delete patches[p];
2539  }
2540 }
2541
2543 {
2544  Array<KnotVector *> kv(3);
2545  NURBSPatchMap p2g(this);
2546
2547  weights.SetSize(GetNDof());
2548  for (int p = 0; p < GetNP(); p++)
2549  {
2550  p2g.SetPatchDofMap(p, kv);
2551  NURBSPatch &Patch = *patches[p];
2552
2553  for (int k = 0; k < kv[2]->GetNCP(); k++)
2554  {
2555  for (int j = 0; j < kv[1]->GetNCP(); j++)
2556  {
2557  for (int i = 0; i < kv[0]->GetNCP(); i++)
2558  {
2559  int l = p2g(i,j,k);
2560  for (int d = 0; d < 3; d++)
2561  {
2562  coords(l*3 + d) = Patch(i,j,k,d)/Patch(i,j,k,3);
2563  }
2564  weights(l) = Patch(i,j,k,3);
2565  }
2566  }
2567  }
2568  delete patches[p];
2569  }
2570 }
2571
2572
2573 #ifdef MFEM_USE_MPI
2575  int *part, const Array<bool> &active_bel)
2576  : gtopo(comm)
2577 {
2578  if (parent->NumOfActiveElems < parent->NumOfElements)
2579  // SetActive (BuildGroups?) and the way the weights are copied
2580  // do not support this case
2581  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
2582  " all elements in the parent must be active!");
2583
2584  patchTopo = parent->patchTopo;
2585  // steal ownership of patchTopo from the 'parent' NURBS extension
2586  if (!parent->own_topo)
2587  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
2588  " parent does not own the patch topology!");
2589  own_topo = 1;
2590  parent->own_topo = 0;
2591
2592  parent->edge_to_knot.Copy(edge_to_knot);
2593
2594  Order = parent->GetOrder();
2595
2596  NumOfKnotVectors = parent->GetNKV();
2597  knotVectors.SetSize(NumOfKnotVectors);
2598  for (int i = 0; i < NumOfKnotVectors; i++)
2599  {
2600  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
2601  }
2602
2603  GenerateOffsets();
2604  CountElements();
2605  CountBdrElements();
2606
2607  // copy 'part' to 'partitioning'
2608  partitioning = new int[GetGNE()];
2609  for (int i = 0; i < GetGNE(); i++)
2610  {
2611  partitioning[i] = part[i];
2612  }
2613  SetActive(partitioning, active_bel);
2614
2617  // GenerateActiveBdrElems(); // done by SetActive for now
2619
2620  Table *serial_elem_dof = parent->GetElementDofTable();
2621  BuildGroups(partitioning, *serial_elem_dof);
2622
2623  weights.SetSize(GetNDof());
2624  // copy weights from parent
2625  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
2626  {
2627  if (activeElem[gel])
2628  {
2629  int ndofs = el_dof->RowSize(lel);
2630  int *ldofs = el_dof->GetRow(lel);
2631  int *gdofs = serial_elem_dof->GetRow(gel);
2632  for (int i = 0; i < ndofs; i++)
2633  {
2634  weights(ldofs[i]) = parent->weights(gdofs[i]);
2635  }
2636  lel++;
2637  }
2638  }
2639 }
2640
2642  ParNURBSExtension *par_parent)
2643  : gtopo(par_parent->gtopo.GetComm())
2644 {
2645  // steal all data from parent
2646  Order = parent->Order;
2647
2648  patchTopo = parent->patchTopo;
2649  own_topo = parent->own_topo;
2650  parent->own_topo = 0;
2651
2652  Swap(edge_to_knot, parent->edge_to_knot);
2653
2655  Swap(knotVectors, parent->knotVectors);
2656
2657  NumOfVertices = parent->NumOfVertices;
2658  NumOfElements = parent->NumOfElements;
2660  NumOfDofs = parent->NumOfDofs;
2661
2662  Swap(v_meshOffsets, parent->v_meshOffsets);
2663  Swap(e_meshOffsets, parent->e_meshOffsets);
2664  Swap(f_meshOffsets, parent->f_meshOffsets);
2665  Swap(p_meshOffsets, parent->p_meshOffsets);
2666
2671
2675  NumOfActiveDofs = parent->NumOfActiveDofs;
2676
2677  Swap(activeVert, parent->activeVert);
2678  Swap(activeElem, parent->activeElem);
2679  Swap(activeBdrElem, parent->activeBdrElem);
2680  Swap(activeDof, parent->activeDof);
2681
2682  el_dof = parent->el_dof;
2683  bel_dof = parent->bel_dof;
2684  parent->el_dof = parent->bel_dof = NULL;
2685
2686  Swap(el_to_patch, parent->el_to_patch);
2687  Swap(bel_to_patch, parent->bel_to_patch);
2688  Swap(el_to_IJK, parent->el_to_IJK);
2689  Swap(bel_to_IJK, parent->bel_to_IJK);
2690
2691  Swap(weights, parent->weights);
2692
2693  delete parent;
2694
2695  partitioning = NULL;
2696
2697 #ifdef MFEM_DEBUG
2698  if (par_parent->partitioning == NULL)
2699  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
2700  " parent ParNURBSExtension has no partitioning!");
2701 #endif
2702
2703  Table *serial_elem_dof = GetGlobalElementDofTable();
2704  BuildGroups(par_parent->partitioning, *serial_elem_dof);
2705  delete serial_elem_dof;
2706 }
2707
2708 Table *ParNURBSExtension::GetGlobalElementDofTable()
2709 {
2710  if (Dimension() == 2)
2711  {
2712  return Get2DGlobalElementDofTable();
2713  }
2714  else
2715  {
2716  return Get3DGlobalElementDofTable();
2717  }
2718 }
2719
2720 Table *ParNURBSExtension::Get2DGlobalElementDofTable()
2721 {
2722  int el = 0;
2723  KnotVector *kv[2];
2724  NURBSPatchMap p2g(this);
2725
2726  Table *gel_dof = new Table(GetGNE(), (Order + 1)*(Order + 1));
2727
2728  for (int p = 0; p < GetNP(); p++)
2729  {
2730  p2g.SetPatchDofMap(p, kv);
2731
2733  for (int j = 0; j < kv[1]->GetNKS(); j++)
2734  {
2735  if (kv[1]->isElement(j))
2736  {
2737  for (int i = 0; i < kv[0]->GetNKS(); i++)
2738  {
2739  if (kv[0]->isElement(i))
2740  {
2741  int *dofs = gel_dof->GetRow(el);
2742  int idx = 0;
2743  for (int jj = 0; jj <= Order; jj++)
2744  {
2745  for (int ii = 0; ii <= Order; ii++)
2746  {
2747  dofs[idx] = p2g(i+ii,j+jj);
2748  idx++;
2749  }
2750  }
2751  el++;
2752  }
2753  }
2754  }
2755  }
2756  }
2757  return gel_dof;
2758 }
2759
2760 Table *ParNURBSExtension::Get3DGlobalElementDofTable()
2761 {
2762  int el = 0;
2763  KnotVector *kv[3];
2764  NURBSPatchMap p2g(this);
2765
2766  Table *gel_dof = new Table(GetGNE(), (Order + 1)*(Order + 1)*(Order + 1));
2767
2768  for (int p = 0; p < GetNP(); p++)
2769  {
2770  p2g.SetPatchDofMap(p, kv);
2771
2773  for (int k = 0; k < kv[2]->GetNKS(); k++)
2774  {
2775  if (kv[2]->isElement(k))
2776  {
2777  for (int j = 0; j < kv[1]->GetNKS(); j++)
2778  {
2779  if (kv[1]->isElement(j))
2780  {
2781  for (int i = 0; i < kv[0]->GetNKS(); i++)
2782  {
2783  if (kv[0]->isElement(i))
2784  {
2785  int *dofs = gel_dof->GetRow(el);
2786  int idx = 0;
2787  for (int kk = 0; kk <= Order; kk++)
2788  {
2789  for (int jj = 0; jj <= Order; jj++)
2790  {
2791  for (int ii = 0; ii <= Order; ii++)
2792  {
2793  dofs[idx] = p2g(i+ii,j+jj,k+kk);
2794  idx++;
2795  }
2796  }
2797  }
2798  el++;
2799  }
2800  }
2801  }
2802  }
2803  }
2804  }
2805  }
2806  return gel_dof;
2807 }
2808
2809 void ParNURBSExtension::SetActive(int *_partitioning,
2810  const Array<bool> &active_bel)
2811 {
2813  activeElem = false;
2814  NumOfActiveElems = 0;
2815  int MyRank = gtopo.MyRank();
2816  for (int i = 0; i < GetGNE(); i++)
2817  if (_partitioning[i] == MyRank)
2818  {
2819  activeElem[i] = true;
2820  NumOfActiveElems++;
2821  }
2822
2823  active_bel.Copy(activeBdrElem);
2824  NumOfActiveBdrElems = 0;
2825  for (int i = 0; i < GetGNBE(); i++)
2826  if (activeBdrElem[i])
2827  {
2829  }
2830 }
2831
2832 void ParNURBSExtension::BuildGroups(int *_partitioning, const Table &elem_dof)
2833 {
2834  Table dof_proc;
2835
2836  ListOfIntegerSets groups;
2837  IntegerSet group;
2838
2839  Transpose(elem_dof, dof_proc); // dof_proc is dof_elem
2840  // convert elements to processors
2841  for (int i = 0; i < dof_proc.Size_of_connections(); i++)
2842  {
2843  dof_proc.GetJ()[i] = _partitioning[dof_proc.GetJ()[i]];
2844  }
2845
2846  // the first group is the local one
2847  int MyRank = gtopo.MyRank();
2848  group.Recreate(1, &MyRank);
2849  groups.Insert(group);
2850
2851  int dof = 0;
2853  for (int d = 0; d < GetNTotalDof(); d++)
2854  if (activeDof[d])
2855  {
2856  group.Recreate(dof_proc.RowSize(d), dof_proc.GetRow(d));
2857  ldof_group[dof] = groups.Insert(group);
2858
2859  dof++;
2860  }
2861
2862  gtopo.Create(groups, 1822);
2863 }
2864 #endif
2865
2866
2867 void NURBSPatchMap::GetPatchKnotVectors(int p, KnotVector *kv[])
2868 {
2869  Ext->patchTopo->GetElementVertices(p, verts);
2870  Ext->patchTopo->GetElementEdges(p, edges, oedge);
2871  if (Ext->Dimension() == 2)
2872  {
2873  kv[0] = Ext->KnotVec(edges[0]);
2874  kv[1] = Ext->KnotVec(edges[1]);
2875  }
2876  else
2877  {
2878  Ext->patchTopo->GetElementFaces(p, faces, oface);
2879
2880  kv[0] = Ext->KnotVec(edges[0]);
2881  kv[1] = Ext->KnotVec(edges[3]);
2882  kv[2] = Ext->KnotVec(edges[8]);
2883  }
2884  opatch = 0;
2885 }
2886
2887 void NURBSPatchMap::GetBdrPatchKnotVectors(int p, KnotVector *kv[], int *okv)
2888 {
2889  Ext->patchTopo->GetBdrElementVertices(p, verts);
2890  Ext->patchTopo->GetBdrElementEdges(p, edges, oedge);
2891  kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]);
2892  if (Ext->Dimension() == 3)
2893  {
2894  faces.SetSize(1);
2895  Ext->patchTopo->GetBdrElementFace(p, &faces[0], &opatch);
2896  kv[1] = Ext->KnotVec(edges[1], oedge[1], &okv[1]);
2897  }
2898  else
2899  {
2900  opatch = oedge[0];
2901  }
2902 }
2903
2905 {
2906  GetPatchKnotVectors(p, kv);
2907
2908  I = kv[0]->GetNE() - 1;
2909  J = kv[1]->GetNE() - 1;
2910
2911  for (int i = 0; i < verts.Size(); i++)
2912  {
2913  verts[i] = Ext->v_meshOffsets[verts[i]];
2914  }
2915
2916  for (int i = 0; i < edges.Size(); i++)
2917  {
2918  edges[i] = Ext->e_meshOffsets[edges[i]];
2919  }
2920
2921  if (Ext->Dimension() == 3)
2922  {
2923  K = kv[2]->GetNE() - 1;
2924
2925  for (int i = 0; i < faces.Size(); i++)
2926  {
2927  faces[i] = Ext->f_meshOffsets[faces[i]];
2928  }
2929  }
2930
2931  pOffset = Ext->p_meshOffsets[p];
2932 }
2933
2935 {
2936  GetPatchKnotVectors(p, kv);
2937
2938  I = kv[0]->GetNCP() - 2;
2939  J = kv[1]->GetNCP() - 2;
2940
2941  for (int i = 0; i < verts.Size(); i++)
2942  {
2943  verts[i] = Ext->v_spaceOffsets[verts[i]];
2944  }
2945
2946  for (int i = 0; i < edges.Size(); i++)
2947  {
2948  edges[i] = Ext->e_spaceOffsets[edges[i]];
2949  }
2950
2951  if (Ext->Dimension() == 3)
2952  {
2953  K = kv[2]->GetNCP() - 2;
2954
2955  for (int i = 0; i < faces.Size(); i++)
2956  {
2957  faces[i] = Ext->f_spaceOffsets[faces[i]];
2958  }
2959  }
2960
2961  pOffset = Ext->p_spaceOffsets[p];
2962 }
2963
2965 {
2966  GetBdrPatchKnotVectors(p, kv, okv);
2967
2968  I = kv[0]->GetNE() - 1;
2969
2970  for (int i = 0; i < verts.Size(); i++)
2971  {
2972  verts[i] = Ext->v_meshOffsets[verts[i]];
2973  }
2974
2975  if (Ext->Dimension() == 2)
2976  {
2977  pOffset = Ext->e_meshOffsets[edges[0]];
2978  }
2979  else
2980  {
2981  J = kv[1]->GetNE() - 1;
2982
2983  for (int i = 0; i < edges.Size(); i++)
2984  {
2985  edges[i] = Ext->e_meshOffsets[edges[i]];
2986  }
2987
2988  pOffset = Ext->f_meshOffsets[faces[0]];
2989  }
2990 }
2991
2992 void NURBSPatchMap::SetBdrPatchDofMap(int p, KnotVector *kv[], int *okv)
2993 {
2994  GetBdrPatchKnotVectors(p, kv, okv);
2995
2996  I = kv[0]->GetNCP() - 2;
2997
2998  for (int i = 0; i < verts.Size(); i++)
2999  {
3000  verts[i] = Ext->v_spaceOffsets[verts[i]];
3001  }
3002
3003  if (Ext->Dimension() == 2)
3004  {
3005  pOffset = Ext->e_spaceOffsets[edges[0]];
3006  }
3007  else
3008  {
3009  J = kv[1]->GetNCP() - 2;
3010
3011  for (int i = 0; i < edges.Size(); i++)
3012  {
3013  edges[i] = Ext->e_spaceOffsets[edges[i]];
3014  }
3015
3016  pOffset = Ext->f_spaceOffsets[faces[0]];
3017  }
3018 }
3019
3020 }
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:2473
void Set3DSolutionVector(Vector &Nodes)
Definition: nurbs.cpp:2542
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:2342
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:1943
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:4084
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:2934
void Generate2DElementDofTable()
Definition: nurbs.cpp:2042
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Definition: nurbs.cpp:1568
Array< int > ldof_group
Definition: nurbs.hpp:350
int GetPatch() const
Definition: fe.hpp:2021
void Generate3DElementDofTable()
Definition: nurbs.cpp:2093
void KnotInsert(int dir, const KnotVector &knot)
Definition: nurbs.cpp:543
void MergeWeights(Mesh *mesh_array[], int num_pieces)
Definition: nurbs.cpp:1543
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:259
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:161
void UniformRefinement()
Definition: nurbs.cpp:2399
void GetElementFaces(int i, Array< int > &, Array< int > &) const
Return the indices and the orientations of all faces of element i.
Definition: mesh.cpp:4281
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:84
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:2434
void Generate3DBdrElementDofTable()
Definition: nurbs.cpp:2216
int NumOfControlPoints
Definition: nurbs.hpp:35
Array< int > e_spaceOffsets
Definition: nurbs.hpp:185
void GenerateBdrElementDofTable()
Definition: nurbs.cpp:2156
void GetElementVertices(int i, Array< int > &dofs) const
Returns the indices of the dofs of element i.
Definition: mesh.hpp:501
const KnotVector * GetKnotVector(int i) const
Definition: nurbs.hpp:300
void Get2DElementTopo(Array< Element * > &elements)
Definition: nurbs.cpp:1848
void GetElementLocalToGlobal(Array< int > &lelem_elem)
Definition: nurbs.cpp:2281
void init(int dim_)
Definition: nurbs.cpp:277
Array< int > edge_to_knot
Definition: nurbs.hpp:173
void GenerateActiveVertices()
Definition: nurbs.cpp:1454
Array< bool > activeElem
Definition: nurbs.hpp:167
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
Definition: nurbs.cpp:2271
Vector knot
Definition: nurbs.hpp:34
int GetElement() const
Definition: fe.hpp:2023
void skip_comment_lines(std::istream &is, const char comment_char)
Definition: mesh.cpp:2249
Array< int > el_to_patch
Definition: nurbs.hpp:191
void SetBdrPatchVertexMap(int p, KnotVector *kv[], int *okv)
Definition: nurbs.cpp:2964
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:1884
KnotVector()
Create KnotVector.
Definition: nurbs.hpp:39
void Get3DBdrElementTopo(Array< Element * > &boundary)
Definition: nurbs.cpp:1974
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:2291
void GetBdrElementTopo(Array< Element * > &boundary)
Definition: nurbs.cpp:1929
void GetElementTopo(Array< Element * > &elements)
Definition: nurbs.cpp:1834
void Generate2DBdrElementDofTable()
Definition: nurbs.cpp:2175
Data type hexahedron element.
Definition: hexahedron.hpp:22
int dim
Definition: ex3.cpp:48
void SetPatch(int p) const
Definition: fe.hpp:2022
void GetPatchKnotVectors(int p, Array< KnotVector * > &kv)
Definition: nurbs.cpp:1667
void GenerateActiveBdrElems()
Definition: nurbs.cpp:1521
double * data
Definition: nurbs.hpp:87
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:140
void LoadPatchTopo(std::istream &input, Array< int > &edge_to_knot)
Definition: mesh.cpp:3582
virtual ~NURBSExtension()
Destroy a NURBSExtension.
Definition: nurbs.cpp:1370
void GenerateOffsets()
Definition: nurbs.cpp:1707
Array< int > activeDof
Definition: nurbs.hpp:169
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:385
static void Get3DRotationMatrix(double n[], double angle, double r, DenseMatrix &T)
Definition: nurbs.cpp:968
void SetKnotsFromPatches()
Definition: nurbs.cpp:2350
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:67
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:270
void SetElement(int e) const
Definition: fe.hpp:2024
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:4319
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:314
Abstract finite element space.
Definition: fespace.hpp:62
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:323
void GetBdrPatchKnotVectors(int p, Array< KnotVector * > &kv)
Definition: nurbs.cpp:1688
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:2515
Table * GetElementDofTable()
Definition: nurbs.hpp:308
Array2D< int > el_to_IJK
Definition: nurbs.hpp:193
void CountBdrElements()
Definition: nurbs.cpp:1814
NURBSExtension * NURBSext
Definition: mesh.hpp:142
ParNURBSExtension(MPI_Comm comm, NURBSExtension *parent, int *partitioning, const Array< bool > &active_bel)
Definition: nurbs.cpp:2574
Array< KnotVector * > & KnotVectors() const
Definition: fe.hpp:2025
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:4062
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:2012
void GetBdrElementVertices(int i, Array< int > &dofs) const
Returns the indices of the dofs of boundary element i.
Definition: mesh.hpp:505
void SetIJK(int *IJK) const
Definition: fe.hpp:2020
void SetBdrPatchDofMap(int p, KnotVector *kv[], int *okv)
Definition: nurbs.cpp:2992
void Difference(const KnotVector &kv, Vector &diff) const
Definition: nurbs.cpp:243
Vector & Weights() const
Definition: fe.hpp:2026
void CalcDShape(Vector &grad, int i, double xi)
Definition: nurbs.cpp:158
void DegreeElevate(int t)
Definition: nurbs.cpp:2391
void Get2DPatchNets(const Vector &Nodes)
Definition: nurbs.cpp:2446
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:2311
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:1429
void CheckBdrPatches()
Definition: nurbs.cpp:1639
void SetPatchVertexMap(int p, KnotVector *kv[])
Definition: nurbs.cpp:2904
void FlipDirection(int dir)
Definition: nurbs.cpp:932
void SetSolutionVector(Vector &Nodes)
Definition: nurbs.cpp:2503
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:2331
void KnotInsert(Array< KnotVector * > &kv)
Definition: nurbs.cpp:2407
void Print(std::ostream &out) const
Definition: nurbs.cpp:1394
int SetLoopDirection(int dir)
Definition: nurbs.cpp:466