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