18#if defined(_MSC_VER) && (_MSC_VER < 1800)
20#define copysign _copysign
53 MFEM_ASSERT(continuity.
Size() == (intervals.
Size() + 1),
54 "Incompatible sizes of continuity and intervals.");
56 const int num_knots =
Order * continuity.
Size() - continuity.
Sum();
59 MFEM_ASSERT(num_knots >= 0,
"Invalid continuity vector for order.");
64 for (
int i = 0; i < continuity.
Size(); ++i)
66 const int multiplicity =
Order - continuity[i];
67 MFEM_ASSERT(multiplicity >= 1 && multiplicity <=
Order+1,
68 "Invalid knot multiplicity for order.");
69 for (
int j = 0; j < multiplicity; ++j)
74 if (i < intervals.
Size()) { accum += intervals[i]; }
79 "Insufficient number of knots to define NURBS.");
82 for (
int i = 0; i <
GetNKS(); ++i)
108 " Parent KnotVector order higher than child");
111 const int nOrder =
Order + t;
114 for (
int i = 0; i <= nOrder; i++)
116 (*newkv)[i] =
knot(0);
118 for (
int i = nOrder + 1; i < newkv->
GetNCP(); i++)
120 (*newkv)[i] =
knot(i - t);
122 for (
int i = 0; i <= nOrder; i++)
134 MFEM_VERIFY(rf > 1,
"Refinement factor must be at least 2.");
140 for (
int i = 0; i <
knot.
Size()-1; i++)
144 for (
int m = 1; m < rf; ++m)
146 newknots(j) = m * h * (
knot(i) +
knot(i+1));
175 if (cf < 2) {
return fine; }
178 MFEM_VERIFY(cne > 0 && cne * cf ==
NumOfElements,
"Invalid coarsening factor");
185 for (
int c=0; c<cne; ++c)
191 if (
knot(i) != kprev)
197 fine[fcnt] =
knot(i);
204 MFEM_VERIFY(fcnt == fine.
Size(),
"");
211 MFEM_VERIFY(rf > 1,
"Refinement factor must be at least 2.");
229 for (
int i = 0; i <
knot.
Size()-1; i++)
254 for (j=0; j<rf-1; ++j)
257 newknots(((rf - 1) * i) + j) = ((1.0 - s0) * k0) + (s0 * k1);
259 s0 += s[(rf*i) + j + 1];
286 for (
int i = 1; i <= ns; i++)
302 MFEM_VERIFY(
GetNE(),
"Elements not counted. Use GetElements().");
312 for (
int e = 0; e <
GetNE(); e++, cnt++)
317 for (
int j = 0; j <samples; j++)
323 for (
int d = 0; d <
Order+1; d++) { os<<
"\t"<<shape[d]; }
326 for (
int d = 0; d <
Order+1; d++) { os<<
"\t"<<shape[d]; }
329 for (
int d = 0; d <
Order+1; d++) { os<<
"\t"<<shape[d]; }
342 int ip = (i >= 0) ? (i + p) : (-1 - i + p);
343 real_t u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), saved, tmp;
344 real_t left[MaxOrder+1], right[MaxOrder+1];
347 for (int j = 1; j <= p; ++j)
349 left[j] = u - knot(ip+1-j);
350 right[j] = knot(ip+j) - u;
352 for (int r = 0; r < j; ++r)
354 tmp = shape(r)/(right[r+1] + left[j-r]);
355 shape(r) = saved + right[r+1]*tmp;
356 saved = left[j-r]*tmp;
362// Routine from "The NURBS book
" - 2nd ed - Piegl and Tiller
363// Algorithm A2.3 p. 72
364void KnotVector::CalcDShape(Vector &grad, int i, real_t xi) const
366 int p = Order, rk, pk;
367 int ip = (i >= 0) ? (i + p) : (-1 - i + p);
368 real_t u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), temp, saved, d;
369 real_t ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1], right[MaxOrder+1];
379 for (int j = 1; j <= p; j++)
381 left[j] = u - knot(ip-j+1);
382 right[j] = knot(ip+j) - u;
384 for (int r = 0; r < j; r++)
386 ndu[j][r] = right[r+1] + left[j-r];
387 temp = ndu[r][j-1]/ndu[j][r];
388 ndu[r][j] = saved + right[r+1]*temp;
389 saved = left[j-r]*temp;
394 for (int r = 0; r <= p; ++r)
401 d = ndu[rk][pk]/ndu[p][rk];
405 d -= ndu[r][pk]/ndu[p][r];
412 grad *= p*(knot(ip+1) - knot(ip));
416 grad *= p*(knot(ip) - knot(ip+1));
420// Routine from "The NURBS book
" - 2nd ed - Piegl and Tiller
421// Algorithm A2.3 p. 72
422void KnotVector::CalcDnShape(Vector &gradn, int n, int i, real_t xi) const
424 int p = Order, rk, pk, j1, j2,r,j,k;
425 int ip = (i >= 0) ? (i + p) : (-1 - i + p);
426 real_t u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip);
427 real_t temp, saved, d;
428 real_t a[2][MaxOrder+1],ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1],
439 for (j = 1; j <= p; j++)
441 left[j] = u - knot(ip-j+1);
442 right[j] = knot(ip+j)- u;
445 for (r = 0; r < j; r++)
447 ndu[j][r] = right[r+1] + left[j-r];
448 temp = ndu[r][j-1]/ndu[j][r];
449 ndu[r][j] = saved + right[r+1]*temp;
450 saved = left[j-r]*temp;
455 for (r = 0; r <= p; r++)
460 for (k = 1; k <= n; k++)
467 a[s2][0] = a[s1][0]/ndu[pk+1][rk];
468 d = a[s2][0]*ndu[rk][pk];
489 for (j = j1; j <= j2; j++)
491 a[s2][j] = (a[s1][j] - a[s1][j-1])/ndu[pk+1][rk+j];
492 d += a[s2][j]*ndu[rk+j][pk];
497 a[s2][k] = - a[s1][k-1]/ndu[pk+1][r];
498 d += a[s2][j]*ndu[rk+j][pk];
509 u = (knot(ip+1) - knot(ip));
513 u = (knot(ip) - knot(ip+1));
517 for (k = 1; k <= n-1; k++) { temp *= (p-k)*u; }
519 for (j = 0; j <= p; j++) { gradn[j] *= temp; }
523void KnotVector::FindMaxima(Array<int> &ks, Vector &xi, Vector &u) const
525 Vector shape(Order+1);
526 Vector maxima(GetNCP());
527 real_t arg1, arg2, arg, max1, max2, max;
529 xi.SetSize(GetNCP());
531 ks.SetSize(GetNCP());
532 for (int j = 0; j <GetNCP(); j++)
535 for (int d = 0; d < Order+1; d++)
541 CalcShape(shape, i, arg1);
545 CalcShape(shape, i, arg2);
548 arg = (arg1 + arg2)/2;
549 CalcShape(shape, i, arg);
552 while ( ( max > max1 ) || (max > max2) )
565 arg = (arg1 + arg2)/2;
566 CalcShape(shape, i, arg);
575 u[j] = getKnotLocation(arg, i+Order);
582// Routine from "The NURBS book
" - 2nd ed - Piegl and Tiller
583// Algorithm A9.1 p. 369
584void KnotVector::FindInterpolant(Array<Vector*> &x)
586 int order = GetOrder();
589 // Find interpolation points
590 Vector xi_args, u_args;
592 FindMaxima(i_args,xi_args, u_args);
594 // Assemble collocation matrix
595 Vector shape(order+1);
596 DenseMatrix A(ncp,ncp);
598 for (int i = 0; i < ncp; i++)
600 CalcShape(shape, i_args[i], xi_args[i]);
601 for (int p = 0; p < order+1; p++)
603 A(i,i_args[i] + p) = shape[p];
610 for (int i= 0; i < x.Size(); i++)
617int KnotVector::findKnotSpan(real_t u) const
621 if (u == knot(NumOfControlPoints+Order))
623 mid = NumOfControlPoints;
628 high = NumOfControlPoints + 1;
629 mid = (low + high)/2;
630 while ( (u < knot(mid-1)) || (u > knot(mid)) )
640 mid = (low + high)/2;
646void KnotVector::Difference(const KnotVector &kv, Vector &diff) const
648 if (Order != kv.GetOrder())
651 " Can not compare
knot vectors with different orders!
");
654 int s = kv.Size() - Size();
657 kv.Difference(*this, diff);
663 if (s == 0) { return; }
667 for (int j = 0; j < kv.Size(); j++)
669 if (abs(knot(i) - kv[j]) < 2 * std::numeric_limits<real_t>::epsilon())
681void NURBSPatch::init(int dim)
683 MFEM_ASSERT(dim > 1, "NURBS patch
dimension (including weight) must be
"
690 ni = kv[0]->GetNCP();
695 data = new real_t[ni*Dim];
698 for (int i = 0; i < ni*Dim; i++)
704 else if (kv.Size() == 2)
706 ni = kv[0]->GetNCP();
707 nj = kv[1]->GetNCP();
708 MFEM_ASSERT(ni > 0 && nj > 0, "Invalid
knot vector dimensions.
");
711 data = new real_t[ni*nj*Dim];
714 for (int i = 0; i < ni*nj*Dim; i++)
720 else if (kv.Size() == 3)
722 ni = kv[0]->GetNCP();
723 nj = kv[1]->GetNCP();
724 nk = kv[2]->GetNCP();
725 MFEM_ASSERT(ni > 0 && nj > 0 && nk > 0,
726 "Invalid
knot vector dimensions.
");
728 data = new real_t[ni*nj*nk*Dim];
731 for (int i = 0; i < ni*nj*nk*Dim; i++)
743NURBSPatch::NURBSPatch(const NURBSPatch &orig)
744 : ni(orig.ni), nj(orig.nj), nk(orig.nk), Dim(orig.Dim),
745 data(NULL), kv(orig.kv.Size()), nd(orig.nd), ls(orig.ls), sd(orig.sd)
747 // Allocate and copy data:
748 const int data_size = Dim*ni*nj*((kv.Size() == 2) ? 1 : nk);
749 data = new real_t[data_size];
750 std::memcpy(data, orig.data, data_size*sizeof(real_t));
752 // Copy the knot vectors:
753 for (int i = 0; i < kv.Size(); i++)
755 kv[i] = new KnotVector(*orig.kv[i]);
759NURBSPatch::NURBSPatch(std::istream &input)
761 int pdim, dim, size = 1;
764 input >> ws >> ident >> pdim; // knotvectors
766 for (int i = 0; i < pdim; i++)
768 kv[i] = new KnotVector(input);
769 size *= kv[i]->GetNCP();
772 input >> ws >> ident >> dim; // dimension
775 input >> ws >> ident; // controlpoints (homogeneous coordinates)
776 if (ident == "controlpoints
" || ident == "controlpoints_homogeneous
")
778 for (int j = 0, i = 0; i < size; i++)
779 for (int d = 0; d <= dim; d++, j++)
784 else // "controlpoints_cartesian
" (Cartesian coordinates with weight)
786 for (int j = 0, i = 0; i < size; i++)
788 for (int d = 0; d <= dim; d++)
792 for (int d = 0; d < dim; d++)
794 data[j+d] *= data[j+dim];
801NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim)
804 kv[0] = new KnotVector(*kv0);
805 kv[1] = new KnotVector(*kv1);
809NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
810 const KnotVector *kv2, int dim)
813 kv[0] = new KnotVector(*kv0);
814 kv[1] = new KnotVector(*kv1);
815 kv[2] = new KnotVector(*kv2);
819NURBSPatch::NURBSPatch(Array<const KnotVector *> &kvs, int dim)
821 kv.SetSize(kvs.Size());
822 for (int i = 0; i < kv.Size(); i++)
824 kv[i] = new KnotVector(*kvs[i]);
829NURBSPatch::NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP)
831 kv.SetSize(parent->kv.Size());
832 for (int i = 0; i < kv.Size(); i++)
835 kv[i] = new KnotVector(*parent->kv[i]);
839 kv[i] = new KnotVector(Order, NCP);
844void NURBSPatch::swap(NURBSPatch *np)
851 for (int i = 0; i < kv.Size(); i++)
853 if (kv[i]) { delete kv[i]; }
870NURBSPatch::~NURBSPatch()
877 for (int i = 0; i < kv.Size(); i++)
879 if (kv[i]) { delete kv[i]; }
883void NURBSPatch::Print(std::ostream &os) const
887 os << "knotvectors\n
" << kv.Size() << '\n';
888 for (int i = 0; i < kv.Size(); i++)
891 size *= kv[i]->GetNCP();
894 os << "\ndimension\n
" << Dim - 1
895 << "\n\ncontrolpoints\n
";
896 for (int j = 0, i = 0; i < size; i++)
899 for (int d = 1; d < Dim; d++)
901 os << ' ' << data[j++];
907int NURBSPatch::SetLoopDirection(int dir)
909 if (nj == -1) // 1D case
920 mfem::err << "NURBSPatch::SetLoopDirection :\n
"
921 " Direction error in 1D patch, dir =
" << dir << '\n';
925 else if (nk == -1) // 2D case
943 mfem::err << "NURBSPatch::SetLoopDirection :\n
"
944 " Direction error in 2D patch, dir =
" << dir << '\n';
973 mfem::err << "NURBSPatch::SetLoopDirection :\n
"
974 " Direction error in 3D patch, dir =
" << dir << '\n';
982void NURBSPatch::UniformRefinement(Array<int> const& rf)
985 for (int dir = 0; dir < kv.Size(); dir++)
989 kv[dir]->Refinement(newknots, rf[dir]);
990 KnotInsert(dir, newknots);
995void NURBSPatch::UniformRefinement(int rf)
997 Array<int> rf_array(kv.Size());
999 UniformRefinement(rf_array);
1002void NURBSPatch::Coarsen(Array<int> const& cf, real_t tol)
1004 for (int dir = 0; dir < kv.Size(); dir++)
1006 if (!kv[dir]->coarse)
1008 const int ne_fine = kv[dir]->GetNE();
1009 KnotRemove(dir, kv[dir]->GetFineKnots(cf[dir]), tol);
1010 kv[dir]->coarse = true;
1011 kv[dir]->GetElements();
1013 const int ne_coarse = kv[dir]->GetNE();
1014 MFEM_VERIFY(ne_fine == cf[dir] * ne_coarse, "");
1015 if (kv[dir]->spacing)
1017 kv[dir]->spacing->SetSize(ne_coarse);
1018 kv[dir]->spacing->ScaleParameters((real_t) cf[dir]);
1024void NURBSPatch::Coarsen(int cf, real_t tol)
1026 Array<int> cf_array(kv.Size());
1028 Coarsen(cf_array, tol);
1031void NURBSPatch::GetCoarseningFactors(Array<int> & f) const
1033 f.SetSize(kv.Size());
1034 for (int dir = 0; dir < kv.Size(); dir++)
1036 f[dir] = kv[dir]->GetCoarseningFactor();
1040void NURBSPatch::KnotInsert(Array<KnotVector *> &newkv)
1042 MFEM_ASSERT(newkv.Size() == kv.Size(), "Invalid input to KnotInsert
");
1043 for (int dir = 0; dir < kv.Size(); dir++)
1045 KnotInsert(dir, *newkv[dir]);
1049void NURBSPatch::KnotInsert(int dir, const KnotVector &newkv)
1051 if (dir >= kv.Size() || dir < 0)
1056 int t = newkv.GetOrder() - kv[dir]->GetOrder();
1060 DegreeElevate(dir, t);
1064 mfem_error("NURBSPatch::KnotInsert : Incorrect order!
");
1068 GetKV(dir)->Difference(newkv, diff);
1069 if (diff.Size() > 0)
1071 KnotInsert(dir, diff);
1075void NURBSPatch::KnotInsert(Array<Vector *> &newkv)
1077 MFEM_ASSERT(newkv.Size() == kv.Size(), "Invalid input to KnotInsert
");
1078 for (int dir = 0; dir < kv.Size(); dir++)
1080 KnotInsert(dir, *newkv[dir]);
1084void NURBSPatch::KnotRemove(Array<Vector *> &rmkv, real_t tol)
1086 for (int dir = 0; dir < kv.Size(); dir++)
1088 KnotRemove(dir, *rmkv[dir], tol);
1092void NURBSPatch::KnotRemove(int dir, const Vector &knot, real_t tol)
1094 // TODO: implement an efficient version of this!
1097 KnotRemove(dir, k, 1, tol);
1101// Algorithm A5.5 from "The NURBS Book
", 2nd ed, Piegl and Tiller, chapter 5.
1102void NURBSPatch::KnotInsert(int dir, const Vector &knot)
1104 if (knot.Size() == 0 ) { return; }
1106 if (dir >= kv.Size() || dir < 0)
1111 NURBSPatch &oldp = *this;
1112 KnotVector &oldkv = *kv[dir];
1114 NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder(),
1115 oldkv.GetNCP() + knot.Size());
1116 NURBSPatch &newp = *newpatch;
1117 KnotVector &newkv = *newp.GetKV(dir);
1119 newkv.spacing = oldkv.spacing;
1121 int size = oldp.SetLoopDirection(dir);
1122 if (size != newp.SetLoopDirection(dir))
1127 int rr = knot.Size() - 1;
1128 int a = oldkv.findKnotSpan(knot(0)) - 1;
1129 int b = oldkv.findKnotSpan(knot(rr)) - 1;
1130 int pl = oldkv.GetOrder();
1131 int ml = oldkv.GetNCP();
1133 for (int j = 0; j <= a; j++)
1135 newkv[j] = oldkv[j];
1137 for (int j = b+pl; j <= ml+pl; j++)
1139 newkv[j+rr+1] = oldkv[j];
1141 for (int k = 0; k <= (a-pl); k++)
1143 for (int ll = 0; ll < size; ll++)
1145 newp.slice(k,ll) = oldp.slice(k,ll);
1148 for (int k = (b-1); k < ml; k++)
1150 for (int ll = 0; ll < size; ll++)
1152 newp.slice(k+rr+1,ll) = oldp.slice(k,ll);
1159 for (int j = rr; j >= 0; j--)
1161 while ( (knot(j) <= oldkv[i]) && (i > a) )
1163 newkv[k] = oldkv[i];
1164 for (int ll = 0; ll < size; ll++)
1166 newp.slice(k-pl-1,ll) = oldp.slice(i-pl-1,ll);
1173 for (int ll = 0; ll < size; ll++)
1175 newp.slice(k-pl-1,ll) = newp.slice(k-pl,ll);
1178 for (int l = 1; l <= pl; l++)
1181 real_t alfa = newkv[k+l] - knot(j);
1182 if (fabs(alfa) == 0.0)
1184 for (int ll = 0; ll < size; ll++)
1186 newp.slice(ind-1,ll) = newp.slice(ind,ll);
1191 alfa = alfa/(newkv[k+l] - oldkv[i-pl+l]);
1192 for (int ll = 0; ll < size; ll++)
1194 newp.slice(ind-1,ll) = alfa*newp.slice(ind-1,ll) +
1195 (1.0-alfa)*newp.slice(ind,ll);
1204 newkv.GetElements();
1209// Algorithm A5.8 from "The NURBS Book
", 2nd ed, Piegl and Tiller, chapter 5.
1210int NURBSPatch::KnotRemove(int dir, real_t knot, int ntimes, real_t tol)
1212 if (dir >= kv.Size() || dir < 0)
1217 NURBSPatch &oldp = *this;
1218 KnotVector &oldkv = *kv[dir];
1220 // Find the index of the last occurrence of the knot.
1222 int multiplicity = 0;
1223 for (int i=0; i<oldkv.Size(); ++i)
1225 if (oldkv[i] == knot)
1232 MFEM_VERIFY(0 < id && id < oldkv.Size() - 1 && ntimes <= multiplicity,
1233 "Only interior knots of sufficient multiplicity may be removed.
");
1235 const int p = oldkv.GetOrder();
1237 NURBSPatch tmpp(this, dir, p, oldkv.GetNCP());
1239 const int size = oldp.SetLoopDirection(dir);
1240 if (size != tmpp.SetLoopDirection(dir))
1246 for (int k = 0; k < oldp.nd; ++k)
1248 for (int ll = 0; ll < size; ll++)
1250 tmpp.slice(k,ll) = oldp.slice(k,ll);
1255 const int s = multiplicity;
1263 Array2D<real_t> temp(last + ntimes + 1, size);
1265 for (int t=0; t<ntimes; ++t)
1267 int off = first - 1; // Difference in index between temp and P.
1269 for (int ll = 0; ll < size; ll++)
1271 temp(0, ll) = oldp.slice(off, ll);
1272 temp(last + 1 - off, ll) = oldp.slice(last + 1, ll);
1276 int jj = last - off;
1280 // Compute new control points for one removal step
1281 const real_t a_i = (knot - oldkv[i]) / (oldkv[i+p+1] - oldkv[i]);
1282 const real_t a_j = (knot - oldkv[j]) / (oldkv[j+p+1] - oldkv[j]);
1284 for (int ll = 0; ll < size; ll++)
1286 temp(ii,ll) = (1.0 / a_i) * oldp.slice(i,ll) -
1287 ((1.0/a_i) - 1.0) * temp(ii - 1, ll);
1289 temp(jj,ll) = (1.0 / (1.0 - a_j)) * (oldp.slice(j,ll) -
1290 (a_j * temp(jj + 1, ll)));
1297 // Check whether knot is removable
1301 for (int ll = 0; ll < size; ll++)
1303 diff[ll] = temp(ii-1, ll) - temp(jj+1, ll);
1308 const real_t a_i = (knot - oldkv[i]) / (oldkv[i+p+1] - oldkv[i]);
1309 for (int ll = 0; ll < size; ll++)
1310 diff[ll] = oldp.slice(i,ll) - (a_i * temp(ii+1, ll))
1311 - ((1.0 - a_i) * temp(ii-1, ll));
1314 const real_t dist = diff.Norml2();
1317 // Removal failed. Return the number of successful removals.
1318 mfem::out << "Knot removal failed after
" << t
1319 << " successful removals
" << endl;
1323 // Note that the new weights may not be positive.
1325 // Save new control points
1331 for (int ll = 0; ll < size; ll++)
1333 tmpp.slice(i,ll) = temp(i - off,ll);
1334 tmpp.slice(j,ll) = temp(j - off,ll);
1342 } // End of loop (t) over ntimes.
1344 const int fout = ((2*r) - s - p) / 2; // First control point out
1348 for (int k=1; k<ntimes; ++k)
1360 NURBSPatch *newpatch = new NURBSPatch(this, dir, p,
1361 oldkv.GetNCP() - ntimes);
1362 NURBSPatch &newp = *newpatch;
1363 if (size != newp.SetLoopDirection(dir))
1368 for (int k = 0; k < fout; ++k)
1370 for (int ll = 0; ll < size; ll++)
1372 newp.slice(k,ll) = oldp.slice(k,ll); // Copy old data
1376 for (int k = i+1; k < oldp.nd; ++k)
1378 for (int ll = 0; ll < size; ll++)
1380 newp.slice(j,ll) = tmpp.slice(k,ll); // Shift
1386 KnotVector &newkv = *newp.GetKV(dir);
1387 MFEM_VERIFY(newkv.Size() == oldkv.Size() - ntimes, "");
1389 newkv.spacing = oldkv.spacing;
1390 newkv.coarse = oldkv.coarse;
1392 for (int k = 0; k < id - ntimes + 1; k++)
1394 newkv[k] = oldkv[k];
1396 for (int k = id + 1; k < oldkv.Size(); k++)
1398 newkv[k - ntimes] = oldkv[k];
1401 newkv.GetElements();
1408void NURBSPatch::DegreeElevate(int t)
1410 for (int dir = 0; dir < kv.Size(); dir++)
1412 DegreeElevate(dir, t);
1416// Routine from "The NURBS book
" - 2nd ed - Piegl and Tiller
1417void NURBSPatch::DegreeElevate(int dir, int t)
1419 if (dir >= kv.Size() || dir < 0)
1424 MFEM_ASSERT(t >= 0, "DegreeElevate cannot decrease the degree.
");
1426 int i, j, k, kj, mpi, mul, mh, kind, cind, first, last;
1427 int r, a, b, oldr, save, s, tr, lbz, rbz, l;
1428 real_t inv, ua, ub, numer, alf, den, bet, gam;
1430 NURBSPatch &oldp = *this;
1431 KnotVector &oldkv = *kv[dir];
1432 oldkv.GetElements();
1434 NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder() + t,
1435 oldkv.GetNCP() + oldkv.GetNE()*t);
1436 NURBSPatch &newp = *newpatch;
1437 KnotVector &newkv = *newp.GetKV(dir);
1439 if (oldkv.spacing) { newkv.spacing = oldkv.spacing; }
1441 int size = oldp.SetLoopDirection(dir);
1442 if (size != newp.SetLoopDirection(dir))
1447 int p = oldkv.GetOrder();
1448 int n = oldkv.GetNCP()-1;
1450 DenseMatrix bezalfs (p+t+1, p+1);
1451 DenseMatrix bpts (p+1, size);
1452 DenseMatrix ebpts (p+t+1, size);
1453 DenseMatrix nextbpts(p-1, size);
1454 Vector alphas (p-1);
1461 Array2D<int> binom(ph+1, ph+1);
1462 for (i = 0; i <= ph; i++)
1464 binom(i,0) = binom(i,i) = 1;
1465 for (j = 1; j < i; j++)
1467 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
1472 bezalfs(ph,p) = 1.0;
1474 for (i = 1; i <= ph2; i++)
1476 inv = 1.0/binom(ph,i);
1478 for (j = max(0,i-t); j <= mpi; j++)
1480 bezalfs(i,j) = inv*binom(p,j)*binom(t,i-j);
1485 for (i = ph2+1; i < ph; i++)
1488 for (j = max(0,i-t); j <= mpi; j++)
1490 bezalfs(i,j) = bezalfs(ph-i,p-j);
1501 for (l = 0; l < size; l++)
1503 newp.slice(0,l) = oldp.slice(0,l);
1505 for (i = 0; i <= ph; i++)
1510 for (i = 0; i <= p; i++)
1512 for (l = 0; l < size; l++)
1514 bpts(i,l) = oldp.slice(i,l);
1521 while (b < m && oldkv[b] == oldkv[b+1]) { b++; }
1529 if (oldr > 0) { lbz = (oldr+2)/2; }
1532 if (r > 0) { rbz = ph-(r+1)/2; }
1538 for (k = p ; k > mul; k--)
1540 alphas[k-mul-1] = numer/(oldkv[a+k]-ua);
1543 for (j = 1; j <= r; j++)
1547 for (k = p; k >= s; k--)
1549 for (l = 0; l < size; l++)
1550 bpts(k,l) = (alphas[k-s]*bpts(k,l) +
1551 (1.0-alphas[k-s])*bpts(k-1,l));
1553 for (l = 0; l < size; l++)
1555 nextbpts(save,l) = bpts(p,l);
1560 for (i = lbz; i <= ph; i++)
1562 for (l = 0; l < size; l++)
1567 for (j = max(0,i-t); j <= mpi; j++)
1569 for (l = 0; l < size; l++)
1571 ebpts(i,l) += bezalfs(i,j)*bpts(j,l);
1581 bet = (ub-newkv[kind-1])/den;
1583 for (tr = 1; tr < oldr; tr++)
1592 alf = (ub-newkv[i])/(ua-newkv[i]);
1593 for (l = 0; l < size; l++)
1595 newp.slice(i,l) = alf*newp.slice(i,l)-(1.0-alf)*newp.slice(i-1,l);
1600 if ((j-tr) <= (kind-ph+oldr))
1602 gam = (ub-newkv[j-tr])/den;
1603 for (l = 0; l < size; l++)
1605 ebpts(kj,l) = gam*ebpts(kj,l) + (1.0-gam)*ebpts(kj+1,l);
1610 for (l = 0; l < size; l++)
1612 ebpts(kj,l) = bet*ebpts(kj,l) + (1.0-bet)*ebpts(kj+1,l);
1627 for (i = 0; i < (ph-oldr); i++)
1633 for (j = lbz; j <= rbz; j++)
1635 for (l = 0; l < size; l++)
1637 newp.slice(cind,l) = ebpts(j,l);
1644 for (j = 0; j <r; j++)
1645 for (l = 0; l < size; l++)
1647 bpts(j,l) = nextbpts(j,l);
1650 for (j = r; j <= p; j++)
1651 for (l = 0; l < size; l++)
1653 bpts(j,l) = oldp.slice(b-p+j,l);
1662 for (i = 0; i <= ph; i++)
1668 newkv.GetElements();
1673void NURBSPatch::FlipDirection(int dir)
1675 int size = SetLoopDirection(dir);
1677 for (int id = 0; id < nd/2; id++)
1678 for (int i = 0; i < size; i++)
1680 Swap<real_t>((*this).slice(id,i), (*this).slice(nd-1-id,i));
1685void NURBSPatch::SwapDirections(int dir1, int dir2)
1687 if (abs(dir1-dir2) == 2)
1690 " directions 0 and 2 are not supported!
");
1693 Array<const KnotVector *> nkv(kv);
1695 Swap<const KnotVector *>(nkv[dir1], nkv[dir2]);
1696 NURBSPatch *newpatch = new NURBSPatch(nkv, Dim);
1698 int size = SetLoopDirection(dir1);
1699 newpatch->SetLoopDirection(dir2);
1701 for (int id = 0; id < nd; id++)
1702 for (int i = 0; i < size; i++)
1704 (*newpatch).slice(id,i) = (*this).slice(id,i);
1710void NURBSPatch::Rotate(real_t angle, real_t n[])
1720 mfem_error("NURBSPatch::Rotate : Specify an angle for
a 3D rotation.
");
1727void NURBSPatch::Get2DRotationMatrix(real_t angle, DenseMatrix &T)
1729 real_t s = sin(angle);
1730 real_t c = cos(angle);
1739void NURBSPatch::Rotate2D(real_t angle)
1747 Vector x(2), y(NULL, 2);
1749 Get2DRotationMatrix(angle, T);
1752 for (int i = 0; i < kv.Size(); i++)
1754 size *= kv[i]->GetNCP();
1757 for (int i = 0; i < size; i++)
1759 y.SetData(data + i*Dim);
1765void NURBSPatch::Get3DRotationMatrix(real_t n[], real_t angle, real_t r,
1769 const real_t l2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
1770 const real_t l = sqrt(l2);
1772 MFEM_ASSERT(l2 > 0.0, "3D rotation axis is undefined
");
1774 if (fabs(angle) == (real_t)(M_PI_2))
1776 s = r*copysign(1., angle);
1780 else if (fabs(angle) == (real_t)(M_PI))
1795 T(0,0) = (n[0]*n[0] + (n[1]*n[1] + n[2]*n[2])*c)/l2;
1796 T(0,1) = -(n[0]*n[1]*c1)/l2 - (n[2]*s)/l;
1797 T(0,2) = -(n[0]*n[2]*c1)/l2 + (n[1]*s)/l;
1798 T(1,0) = -(n[0]*n[1]*c1)/l2 + (n[2]*s)/l;
1799 T(1,1) = (n[1]*n[1] + (n[0]*n[0] + n[2]*n[2])*c)/l2;
1800 T(1,2) = -(n[1]*n[2]*c1)/l2 - (n[0]*s)/l;
1801 T(2,0) = -(n[0]*n[2]*c1)/l2 - (n[1]*s)/l;
1802 T(2,1) = -(n[1]*n[2]*c1)/l2 + (n[0]*s)/l;
1803 T(2,2) = (n[2]*n[2] + (n[0]*n[0] + n[1]*n[1])*c)/l2;
1806void NURBSPatch::Rotate3D(real_t n[], real_t angle)
1814 Vector x(3), y(NULL, 3);
1816 Get3DRotationMatrix(n, angle, 1., T);
1819 for (int i = 0; i < kv.Size(); i++)
1821 size *= kv[i]->GetNCP();
1824 for (int i = 0; i < size; i++)
1826 y.SetData(data + i*Dim);
1832int NURBSPatch::MakeUniformDegree(int degree)
1838 for (int dir = 0; dir < kv.Size(); dir++)
1840 maxd = std::max(maxd, kv[dir]->GetOrder());
1844 for (int dir = 0; dir < kv.Size(); dir++)
1846 if (maxd > kv[dir]->GetOrder())
1848 DegreeElevate(dir, maxd - kv[dir]->GetOrder());
1855NURBSPatch *Interpolate(NURBSPatch &p1, NURBSPatch &p2)
1857 if (p1.kv.Size() != p2.kv.Size() || p1.Dim != p2.Dim)
1862 int size = 1, dim = p1.Dim;
1863 Array<const KnotVector *> kv(p1.kv.Size() + 1);
1865 for (int i = 0; i < p1.kv.Size(); i++)
1867 if (p1.kv[i]->GetOrder() < p2.kv[i]->GetOrder())
1869 p1.KnotInsert(i, *p2.kv[i]);
1870 p2.KnotInsert(i, *p1.kv[i]);
1874 p2.KnotInsert(i, *p1.kv[i]);
1875 p1.KnotInsert(i, *p2.kv[i]);
1878 size *= kv[i]->GetNCP();
1881 KnotVector &nkv = *(new KnotVector(1, 2));
1882 nkv[0] = nkv[1] = 0.0;
1883 nkv[2] = nkv[3] = 1.0;
1887 NURBSPatch *patch = new NURBSPatch(kv, dim);
1890 for (int i = 0; i < size; i++)
1892 for (int d = 0; d < dim; d++)
1894 patch->data[i*dim+d] = p1.data[i*dim+d];
1895 patch->data[(i+size)*dim+d] = p2.data[i*dim+d];
1902NURBSPatch *Revolve3D(NURBSPatch &patch, real_t n[], real_t ang, int times)
1910 Array<const KnotVector *> nkv(patch.kv.Size() + 1);
1912 for (int i = 0; i < patch.kv.Size(); i++)
1914 nkv[i] = patch.kv[i];
1915 size *= nkv[i]->GetNCP();
1918 KnotVector &lkv = *(new KnotVector(2, ns));
1920 lkv[0] = lkv[1] = lkv[2] = 0.0;
1921 for (int i = 1; i < times; i++)
1923 lkv[2*i+1] = lkv[2*i+2] = i;
1925 lkv[ns] = lkv[ns+1] = lkv[ns+2] = times;
1927 NURBSPatch *newpatch = new NURBSPatch(nkv, 4);
1930 DenseMatrix T(3), T2(3);
1931 Vector u(NULL, 3), v(NULL, 3);
1933 NURBSPatch::Get3DRotationMatrix(n, ang, 1., T);
1934 real_t c = cos(ang/2);
1935 NURBSPatch::Get3DRotationMatrix(n, ang/2, 1./c, T2);
1938 real_t *op = patch.data, *np;
1939 for (int i = 0; i < size; i++)
1941 np = newpatch->data + 4*i;
1942 for (int j = 0; j < 4; j++)
1946 for (int j = 0; j < times; j++)
1949 v.SetData(np += 4*size);
1952 v.SetData(np += 4*size);
1962void NURBSPatch::SetKnotVectorsCoarse(bool c)
1964 for (int i=0; i<kv.Size(); ++i)
1970NURBSExtension::NURBSExtension(const NURBSExtension &orig)
1971 : mOrder(orig.mOrder), mOrders(orig.mOrders),
1972 NumOfKnotVectors(orig.NumOfKnotVectors),
1973 NumOfVertices(orig.NumOfVertices),
1974 NumOfElements(orig.NumOfElements),
1975 NumOfBdrElements(orig.NumOfBdrElements),
1976 NumOfDofs(orig.NumOfDofs),
1977 NumOfActiveVertices(orig.NumOfActiveVertices),
1978 NumOfActiveElems(orig.NumOfActiveElems),
1979 NumOfActiveBdrElems(orig.NumOfActiveBdrElems),
1980 NumOfActiveDofs(orig.NumOfActiveDofs),
1981 activeVert(orig.activeVert),
1982 activeElem(orig.activeElem),
1983 activeBdrElem(orig.activeBdrElem),
1984 activeDof(orig.activeDof),
1985 patchTopo(new Mesh(*orig.patchTopo)),
1987 edge_to_knot(orig.edge_to_knot),
1988 knotVectors(orig.knotVectors.Size()), // knotVectors are copied in the body
1989 knotVectorsCompr(orig.knotVectorsCompr.Size()),
1990 weights(orig.weights),
1991 d_to_d(orig.d_to_d),
1992 master(orig.master),
1994 v_meshOffsets(orig.v_meshOffsets),
1995 e_meshOffsets(orig.e_meshOffsets),
1996 f_meshOffsets(orig.f_meshOffsets),
1997 p_meshOffsets(orig.p_meshOffsets),
1998 v_spaceOffsets(orig.v_spaceOffsets),
1999 e_spaceOffsets(orig.e_spaceOffsets),
2000 f_spaceOffsets(orig.f_spaceOffsets),
2001 p_spaceOffsets(orig.p_spaceOffsets),
2002 el_dof(orig.el_dof ? new Table(*orig.el_dof) : NULL),
2003 bel_dof(orig.bel_dof ? new Table(*orig.bel_dof) : NULL),
2004 el_to_patch(orig.el_to_patch),
2005 bel_to_patch(orig.bel_to_patch),
2006 el_to_IJK(orig.el_to_IJK),
2007 bel_to_IJK(orig.bel_to_IJK),
2008 patches(orig.patches.Size()) // patches are copied in the body
2010 // Copy the knot vectors:
2011 for (int i = 0; i < knotVectors.Size(); i++)
2013 knotVectors[i] = new KnotVector(*orig.knotVectors[i]);
2015 CreateComprehensiveKV();
2017 // Copy the patches:
2018 for (int p = 0; p < patches.Size(); p++)
2020 patches[p] = new NURBSPatch(*orig.patches[p]);
2024NURBSExtension::NURBSExtension(std::istream &input, bool spacing)
2027 patchTopo = new Mesh;
2028 patchTopo->LoadPatchTopo(input, edge_to_knot);
2032 // CheckBdrPatches();
2034 skip_comment_lines(input, '#');
2036 // Read knotvectors or patches
2038 input >> ws >> ident; // 'knotvectors' or 'patches'
2039 if (ident == "knotvectors
")
2041 input >> NumOfKnotVectors;
2042 knotVectors.SetSize(NumOfKnotVectors);
2043 for (int i = 0; i < NumOfKnotVectors; i++)
2045 knotVectors[i] = new KnotVector(input);
2048 if (spacing) // Read spacing formulas for knotvectors
2050 input >> ws >> ident; // 'spacing'
2051 MFEM_VERIFY(ident == "spacing",
2052 "Spacing formula section missing from NURBS mesh file
");
2054 input >> numSpacing;
2055 for (int j = 0; j < numSpacing; j++)
2057 int ki, spacingType, numIntParam, numRealParam;
2058 input >> ki >> spacingType >> numIntParam >> numRealParam;
2060 MFEM_VERIFY(0 <= ki && ki < NumOfKnotVectors,
2061 "Invalid knotvector
index");
2062 MFEM_VERIFY(numIntParam >= 0 && numRealParam >= 0,
2063 "Invalid number of parameters in
KnotVector");
2065 Array<int> ipar(numIntParam);
2066 Vector dpar(numRealParam);
2068 for (int i=0; i<numIntParam; ++i)
2073 for (int i=0; i<numRealParam; ++i)
2078 const SpacingType s = (SpacingType) spacingType;
2079 knotVectors[ki]->spacing = GetSpacingFunction(s, ipar, dpar);
2083 else if (ident == "patches
")
2085 patches.SetSize(GetNP());
2086 for (int p = 0; p < patches.Size(); p++)
2088 skip_comment_lines(input, '#');
2089 patches[p] = new NURBSPatch(input);
2092 NumOfKnotVectors = 0;
2093 for (int i = 0; i < patchTopo->GetNEdges(); i++)
2094 if (NumOfKnotVectors < KnotInd(i))
2096 NumOfKnotVectors = KnotInd(i);
2099 knotVectors.SetSize(NumOfKnotVectors);
2102 Array<int> edges, oedge;
2103 for (int p = 0; p < patches.Size(); p++)
2105 if (Dimension() == 1)
2107 if (knotVectors[KnotInd(p)] == NULL)
2109 knotVectors[KnotInd(p)] =
2110 new KnotVector(*patches[p]->GetKV(0));
2113 else if (Dimension() == 2)
2115 patchTopo->GetElementEdges(p, edges, oedge);
2116 if (knotVectors[KnotInd(edges[0])] == NULL)
2118 knotVectors[KnotInd(edges[0])] =
2119 new KnotVector(*patches[p]->GetKV(0));
2121 if (knotVectors[KnotInd(edges[1])] == NULL)
2123 knotVectors[KnotInd(edges[1])] =
2124 new KnotVector(*patches[p]->GetKV(1));
2127 else if (Dimension() == 3)
2129 patchTopo->GetElementEdges(p, edges, oedge);
2130 if (knotVectors[KnotInd(edges[0])] == NULL)
2132 knotVectors[KnotInd(edges[0])] =
2133 new KnotVector(*patches[p]->GetKV(0));
2135 if (knotVectors[KnotInd(edges[3])] == NULL)
2137 knotVectors[KnotInd(edges[3])] =
2138 new KnotVector(*patches[p]->GetKV(1));
2140 if (knotVectors[KnotInd(edges[8])] == NULL)
2142 knotVectors[KnotInd(edges[8])] =
2143 new KnotVector(*patches[p]->GetKV(2));
2150 MFEM_ABORT("invalid section:
" << ident);
2153 CreateComprehensiveKV();
2155 SetOrdersFromKnotVectors();
2160 // NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs
2162 skip_comment_lines(input, '#');
2164 // Check for a list of mesh elements
2165 if (patches.Size() == 0)
2167 input >> ws >> ident;
2169 if (patches.Size() == 0 && ident == "mesh_elements
")
2171 input >> NumOfActiveElems;
2172 activeElem.SetSize(GetGNE());
2175 for (int i = 0; i < NumOfActiveElems; i++)
2178 activeElem[glob_elem] = true;
2181 skip_comment_lines(input, '#');
2182 input >> ws >> ident;
2186 NumOfActiveElems = NumOfElements;
2187 activeElem.SetSize(NumOfElements);
2191 GenerateActiveVertices();
2193 GenerateElementDofTable();
2194 GenerateActiveBdrElems();
2195 GenerateBdrElementDofTable();
2198 if (ident == "periodic
")
2203 skip_comment_lines(input, '#');
2204 input >> ws >> ident;
2207 if (patches.Size() == 0)
2210 if (ident == "weights
")
2212 weights.Load(input, GetNDof());
2214 else // e.g. ident = "unitweights
" or "autoweights
"
2216 weights.SetSize(GetNDof());
2222 ConnectBoundaries();
2225NURBSExtension::NURBSExtension(NURBSExtension *parent, int newOrder)
2227 patchTopo = parent->patchTopo;
2230 parent->edge_to_knot.Copy(edge_to_knot);
2232 NumOfKnotVectors = parent->GetNKV();
2233 knotVectors.SetSize(NumOfKnotVectors);
2234 knotVectorsCompr.SetSize(parent->GetNP()*parent->Dimension());
2235 const Array<int> &pOrders = parent->GetOrders();
2236 for (int i = 0; i < NumOfKnotVectors; i++)
2238 if (newOrder > pOrders[i])
2241 parent->GetKnotVector(i)->DegreeElevate(newOrder - pOrders[i]);
2245 knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
2248 CreateComprehensiveKV();
2250 // copy some data from parent
2251 NumOfElements = parent->NumOfElements;
2252 NumOfBdrElements = parent->NumOfBdrElements;
2254 SetOrdersFromKnotVectors();
2256 GenerateOffsets(); // dof offsets will be different from parent
2258 NumOfActiveVertices = parent->NumOfActiveVertices;
2259 NumOfActiveElems = parent->NumOfActiveElems;
2260 NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
2261 parent->activeVert.Copy(activeVert);
2263 parent->activeElem.Copy(activeElem);
2264 parent->activeBdrElem.Copy(activeBdrElem);
2266 GenerateElementDofTable();
2267 GenerateBdrElementDofTable();
2269 weights.SetSize(GetNDof());
2273 parent->master.Copy(master);
2274 parent->slave.Copy(slave);
2275 ConnectBoundaries();
2278NURBSExtension::NURBSExtension(NURBSExtension *parent,
2279 const Array<int> &newOrders, Mode mode)
2282 newOrders.Copy(mOrders);
2283 SetOrderFromOrders();
2285 patchTopo = parent->patchTopo;
2288 parent->edge_to_knot.Copy(edge_to_knot);
2290 NumOfKnotVectors = parent->GetNKV();
2291 MFEM_VERIFY(mOrders.Size() == NumOfKnotVectors, "invalid newOrders array
");
2292 knotVectors.SetSize(NumOfKnotVectors);
2293 const Array<int> &pOrders = parent->GetOrders();
2295 for (int i = 0; i < NumOfKnotVectors; i++)
2297 if (mOrders[i] > pOrders[i])
2300 parent->GetKnotVector(i)->DegreeElevate(mOrders[i] - pOrders[i]);
2304 knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
2307 CreateComprehensiveKV();
2309 // copy some data from parent
2310 NumOfElements = parent->NumOfElements;
2311 NumOfBdrElements = parent->NumOfBdrElements;
2313 GenerateOffsets(); // dof offsets will be different from parent
2315 NumOfActiveVertices = parent->NumOfActiveVertices;
2316 NumOfActiveElems = parent->NumOfActiveElems;
2317 NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
2318 parent->activeVert.Copy(activeVert);
2320 parent->activeElem.Copy(activeElem);
2321 parent->activeBdrElem.Copy(activeBdrElem);
2323 GenerateElementDofTable();
2324 GenerateBdrElementDofTable();
2326 weights.SetSize(GetNDof());
2329 parent->master.Copy(master);
2330 parent->slave.Copy(slave);
2331 ConnectBoundaries();
2334NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
2336 NURBSExtension *parent = mesh_array[0]->NURBSext;
2338 if (!parent->own_topo)
2341 " parent does not own the patch topology!
");
2343 patchTopo = parent->patchTopo;
2345 parent->own_topo = false;
2347 parent->edge_to_knot.Copy(edge_to_knot);
2349 parent->GetOrders().Copy(mOrders);
2350 mOrder = parent->GetOrder();
2352 NumOfKnotVectors = parent->GetNKV();
2353 knotVectors.SetSize(NumOfKnotVectors);
2354 for (int i = 0; i < NumOfKnotVectors; i++)
2356 knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
2358 CreateComprehensiveKV();
2364 // assuming the meshes define a partitioning of all the elements
2365 NumOfActiveElems = NumOfElements;
2366 activeElem.SetSize(NumOfElements);
2369 GenerateActiveVertices();
2371 GenerateElementDofTable();
2372 GenerateActiveBdrElems();
2373 GenerateBdrElementDofTable();
2375 weights.SetSize(GetNDof());
2376 MergeWeights(mesh_array, num_pieces);
2379NURBSExtension::NURBSExtension(const Mesh *patch_topology,
2380 const Array<const NURBSPatch*> p)
2382 patchTopo = new Mesh( *patch_topology );
2383 patchTopo->GetEdgeVertexTable();
2385 patches.Reserve(p.Size());
2389 edge_to_knot.SetSize(patch_topology->GetNEdges());
2390 NumOfKnotVectors = 0;
2392 for (int ielem = 0; ielem < patch_topology->GetNE(); ++ielem)
2394 patches.Append(new NURBSPatch(*p[ielem]));
2395 NURBSPatch& patch = *patches[ielem];
2396 int num_patch_elems = 1;
2397 for (int ikv = 0; ikv < patch.GetNKV(); ++ikv)
2399 kvs[ikv] = knotVectors.Size();
2400 knotVectors.Append(new KnotVector(*patch.GetKV(ikv)));
2401 num_patch_elems *= patch.GetKV(ikv)->GetNE();
2404 NumOfElements += num_patch_elems;
2405 patch_topology->GetElementEdges(ielem, edges, oedges);
2406 for (int iedge = 0; iedge < edges.Size(); ++iedge)
2412 edge_to_knot[edges[iedge]] = kvs[1];
2416 edge_to_knot[edges[iedge]] = kvs[0];
2421 edge_to_knot[edges[iedge]] = kvs[2];
2428 NumOfActiveElems = NumOfElements;
2429 activeElem.SetSize(NumOfElements);
2432 CreateComprehensiveKV();
2433 SetOrdersFromKnotVectors();
2435 GenerateActiveVertices();
2437 GenerateElementDofTable();
2438 GenerateActiveBdrElems();
2439 GenerateBdrElementDofTable();
2441 weights.SetSize(GetNDof());
2446NURBSExtension::~NURBSExtension()
2448 if (bel_dof) { delete bel_dof; }
2449 if (el_dof) { delete el_dof; }
2451 for (int i = 0; i < knotVectors.Size(); i++)
2453 delete knotVectors[i];
2456 for (int i = 0; i < knotVectorsCompr.Size(); i++)
2458 delete knotVectorsCompr[i];
2461 for (int i = 0; i < patches.Size(); i++)
2472void NURBSExtension::Print(std::ostream &os, const std::string &comments) const
2474 Array<int> kvSpacing;
2475 if (patches.Size() == 0)
2477 for (int i = 0; i < NumOfKnotVectors; i++)
2479 if (knotVectors[i]->spacing) { kvSpacing.Append(i); }
2483 const int version = kvSpacing.Size() > 0 ? 11 : 10; // v1.0 or v1.1
2484 patchTopo->PrintTopo(os, edge_to_knot, version, comments);
2485 if (patches.Size() == 0)
2487 os << "\nknotvectors\n
" << NumOfKnotVectors << '\n';
2488 for (int i = 0; i < NumOfKnotVectors; i++)
2490 knotVectors[i]->Print(os);
2493 if (kvSpacing.Size() > 0)
2495 os << "\nspacing\n
" << kvSpacing.Size() << '\n';
2496 for (auto kv : kvSpacing)
2499 knotVectors[kv]->spacing->Print(os);
2503 if (NumOfActiveElems < NumOfElements)
2505 os << "\nmesh_elements\n
" << NumOfActiveElems << '\n';
2506 for (int i = 0; i < NumOfElements; i++)
2513 os << "\nweights\n
";
2514 weights.Print(os, 1);
2518 os << "\npatches\n
";
2519 for (int p = 0; p < patches.Size(); p++)
2521 os << "\n# patch
" << p << "\n\n
";
2522 patches[p]->Print(os);
2527void NURBSExtension::PrintCharacteristics(std::ostream &os) const
2530 "NURBS
Mesh entity sizes:\n
"
2531 "Dimension =
" << Dimension() << "\n
"
2533 Array<int> unique_orders(mOrders);
2534 unique_orders.Sort();
2535 unique_orders.Unique();
2536 unique_orders.Print(os, unique_orders.Size());
2538 "NumOfKnotVectors =
" << GetNKV() << "\n
"
2539 "NumOfPatches =
" << GetNP() << "\n
"
2540 "NumOfBdrPatches =
" << GetNBP() << "\n
"
2541 "NumOfVertices =
" << GetGNV() << "\n
"
2543 "NumOfBdrElements =
" << GetGNBE() << "\n
"
2544 "NumOfDofs =
" << GetNTotalDof() << "\n
"
2545 "NumOfActiveVertices =
" << GetNV() << "\n
"
2546 "NumOfActiveElems =
" << GetNE() << "\n
"
2547 "NumOfActiveBdrElems =
" << GetNBE() << "\n
"
2548 "NumOfActiveDofs =
" << GetNDof() << '\n';
2549 for (int i = 0; i < NumOfKnotVectors; i++)
2551 os << ' ' << i + 1 << ")
";
2552 knotVectors[i]->Print(os);
2557void NURBSExtension::PrintFunctions(const char *basename, int samples) const
2560 for (int i = 0; i < NumOfKnotVectors; i++)
2562 std::ostringstream filename;
2563 filename << basename << "_
" << i << ".dat
";
2564 os.open(filename.str().c_str());
2565 knotVectors[i]->PrintFunctions(os,samples);
2570void NURBSExtension::InitDofMap()
2577void NURBSExtension::ConnectBoundaries(Array<int> &bnds0, Array<int> &bnds1)
2581 ConnectBoundaries();
2584void NURBSExtension::ConnectBoundaries()
2586 if (master.Size() != slave.Size())
2588 mfem_error("NURBSExtension::ConnectBoundaries() boundary lists not of equal size
");
2590 if (master.Size() == 0 ) { return; }
2592 // Initialize d_to_d
2593 d_to_d.SetSize(NumOfDofs);
2594 for (int i = 0; i < NumOfDofs; i++) { d_to_d[i] = i; }
2597 for (int i = 0; i < master.Size(); i++)
2599 int bnd0 = -1, bnd1 = -1;
2600 for (int b = 0; b < GetNBP(); b++)
2602 if (master[i] == patchTopo->GetBdrAttribute(b)) { bnd0 = b; }
2603 if (slave[i]== patchTopo->GetBdrAttribute(b)) { bnd1 = b; }
2605 MFEM_VERIFY(bnd0 != -1,"Bdr 0 not found
");
2606 MFEM_VERIFY(bnd1 != -1,"Bdr 1 not found
");
2608 if (Dimension() == 1)
2610 ConnectBoundaries1D(bnd0, bnd1);
2612 else if (Dimension() == 2)
2614 ConnectBoundaries2D(bnd0, bnd1);
2618 ConnectBoundaries3D(bnd0, bnd1);
2623 Array<int> tmp(d_to_d.Size()+1);
2626 for (int i = 0; i < d_to_d.Size(); i++)
2632 for (int i = 0; i < tmp.Size(); i++)
2634 if (tmp[i] == 1) { tmp[i] = cnt++; }
2638 for (int i = 0; i < d_to_d.Size(); i++)
2640 d_to_d[i] = tmp[d_to_d[i]];
2644 if (el_dof) { delete el_dof; }
2645 if (bel_dof) { delete bel_dof; }
2646 GenerateElementDofTable();
2647 GenerateBdrElementDofTable();
2650void NURBSExtension::ConnectBoundaries1D(int bnd0, int bnd1)
2652 NURBSPatchMap p2g0(this);
2653 NURBSPatchMap p2g1(this);
2655 int okv0[1],okv1[1];
2656 const KnotVector *kv0[1],*kv1[1];
2658 p2g0.SetBdrPatchDofMap(bnd0, kv0, okv0);
2659 p2g1.SetBdrPatchDofMap(bnd1, kv1, okv1);
2661 d_to_d[p2g0(0)] = d_to_d[p2g1(0)];
2664void NURBSExtension::ConnectBoundaries2D(int bnd0, int bnd1)
2666 NURBSPatchMap p2g0(this);
2667 NURBSPatchMap p2g1(this);
2669 int okv0[1],okv1[1];
2670 const KnotVector *kv0[1],*kv1[1];
2672 p2g0.SetBdrPatchDofMap(bnd0, kv0, okv0);
2673 p2g1.SetBdrPatchDofMap(bnd1, kv1, okv1);
2676 int nks0 = kv0[0]->GetNKS();
2679 bool compatible = true;
2680 if (p2g0.nx() != p2g1.nx()) { compatible = false; }
2681 if (kv0[0]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
2682 if (kv0[0]->GetOrder() != kv1[0]->GetOrder()) { compatible = false; }
2686 mfem::out<<p2g0.nx()<<" "<<p2g1.nx()<<endl;
2687 mfem::out<<kv0[0]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
2688 mfem::out<<kv0[0]->GetOrder()<<" "<<kv1[0]->GetOrder()<<endl;
2689 mfem_error("NURBS boundaries not compatible
");
2693 for (int i = 0; i < nks0; i++)
2695 if (kv0[0]->isElement(i))
2697 if (!kv1[0]->isElement(i)) { mfem_error("isElement does not match
"); }
2698 for (int ii = 0; ii <= kv0[0]->GetOrder(); ii++)
2700 int ii0 = (okv0[0] >= 0) ? (i+ii) : (nx-i-ii);
2701 int ii1 = (okv1[0] >= 0) ? (i+ii) : (nx-i-ii);
2703 d_to_d[p2g0(ii0)] = d_to_d[p2g1(ii1)];
2710void NURBSExtension::ConnectBoundaries3D(int bnd0, int bnd1)
2712 NURBSPatchMap p2g0(this);
2713 NURBSPatchMap p2g1(this);
2715 int okv0[2],okv1[2];
2716 const KnotVector *kv0[2],*kv1[2];
2718 p2g0.SetBdrPatchDofMap(bnd0, kv0, okv0);
2719 p2g1.SetBdrPatchDofMap(bnd1, kv1, okv1);
2724 int nks0 = kv0[0]->GetNKS();
2725 int nks1 = kv0[1]->GetNKS();
2728 bool compatible = true;
2729 if (p2g0.nx() != p2g1.nx()) { compatible = false; }
2730 if (p2g0.ny() != p2g1.ny()) { compatible = false; }
2732 if (kv0[0]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
2733 if (kv0[1]->GetNKS() != kv1[1]->GetNKS()) { compatible = false; }
2735 if (kv0[0]->GetOrder() != kv1[0]->GetOrder()) { compatible = false; }
2736 if (kv0[1]->GetOrder() != kv1[1]->GetOrder()) { compatible = false; }
2740 mfem::out<<p2g0.nx()<<" "<<p2g1.nx()<<endl;
2741 mfem::out<<p2g0.ny()<<" "<<p2g1.ny()<<endl;
2743 mfem::out<<kv0[0]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
2744 mfem::out<<kv0[1]->GetNKS()<<" "<<kv1[1]->GetNKS()<<endl;
2746 mfem::out<<kv0[0]->GetOrder()<<" "<<kv1[0]->GetOrder()<<endl;
2747 mfem::out<<kv0[1]->GetOrder()<<" "<<kv1[1]->GetOrder()<<endl;
2748 mfem_error("NURBS boundaries not compatible
");
2752 for (int j = 0; j < nks1; j++)
2754 if (kv0[1]->isElement(j))
2756 if (!kv1[1]->isElement(j)) { mfem_error("isElement does not match #1
"); }
2757 for (int i = 0; i < nks0; i++)
2759 if (kv0[0]->isElement(i))
2761 if (!kv1[0]->isElement(i)) { mfem_error("isElement does not match #0
"); }
2762 for (int jj = 0; jj <= kv0[1]->GetOrder(); jj++)
2764 int jj0 = (okv0[1] >= 0) ? (j+jj) : (ny-j-jj);
2765 int jj1 = (okv1[1] >= 0) ? (j+jj) : (ny-j-jj);
2767 for (int ii = 0; ii <= kv0[0]->GetOrder(); ii++)
2769 int ii0 = (okv0[0] >= 0) ? (i+ii) : (nx-i-ii);
2770 int ii1 = (okv1[0] >= 0) ? (i+ii) : (nx-i-ii);
2772 d_to_d[p2g0(ii0,jj0)] = d_to_d[p2g1(ii1,jj1)];
2781void NURBSExtension::GenerateActiveVertices()
2783 int vert[8], nv, g_el, nx, ny, nz, dim = Dimension();
2785 NURBSPatchMap p2g(this);
2786 const KnotVector *kv[3];
2789 activeVert.SetSize(GetGNV());
2791 for (int p = 0; p < GetNP(); p++)
2793 p2g.SetPatchVertexMap(p, kv);
2796 ny = (dim >= 2) ? p2g.ny() : 1;
2797 nz = (dim == 3) ? p2g.nz() : 1;
2799 for (int k = 0; k < nz; k++)
2801 for (int j = 0; j < ny; j++)
2803 for (int i = 0; i < nx; i++)
2805 if (activeElem[g_el])
2815 vert[0] = p2g(i, j );
2816 vert[1] = p2g(i+1,j );
2817 vert[2] = p2g(i+1,j+1);
2818 vert[3] = p2g(i, j+1);
2823 vert[0] = p2g(i, j, k);
2824 vert[1] = p2g(i+1,j, k);
2825 vert[2] = p2g(i+1,j+1,k);
2826 vert[3] = p2g(i, j+1,k);
2828 vert[4] = p2g(i, j, k+1);
2829 vert[5] = p2g(i+1,j, k+1);
2830 vert[6] = p2g(i+1,j+1,k+1);
2831 vert[7] = p2g(i, j+1,k+1);
2835 for (int v = 0; v < nv; v++)
2837 activeVert[vert[v]] = 1;
2846 NumOfActiveVertices = 0;
2847 for (int i = 0; i < GetGNV(); i++)
2848 if (activeVert[i] == 1)
2850 activeVert[i] = NumOfActiveVertices++;
2854void NURBSExtension::GenerateActiveBdrElems()
2856 int dim = Dimension();
2857 Array<KnotVector *> kv(dim);
2859 activeBdrElem.SetSize(GetGNBE());
2860 if (GetGNE() == GetNE())
2862 activeBdrElem = true;
2863 NumOfActiveBdrElems = GetGNBE();
2866 activeBdrElem = false;
2867 NumOfActiveBdrElems = 0;
2868 // the mesh will generate the actual boundary including boundary
2869 // elements that are not on boundary patches. we use this for
2870 // visualization of processor boundaries
2872 // TODO: generate actual boundary?
2876void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces)
2878 Array<int> lelem_elem;
2880 for (int i = 0; i < num_pieces; i++)
2882 NURBSExtension *lext = mesh_array[i]->NURBSext;
2884 lext->GetElementLocalToGlobal(lelem_elem);
2886 for (int lel = 0; lel < lext->GetNE(); lel++)
2888 int gel = lelem_elem[lel];
2890 int nd = el_dof->RowSize(gel);
2891 int *gdofs = el_dof->GetRow(gel);
2892 int *ldofs = lext->el_dof->GetRow(lel);
2893 for (int j = 0; j < nd; j++)
2895 weights(gdofs[j]) = lext->weights(ldofs[j]);
2901void NURBSExtension::MergeGridFunctions(
2902 GridFunction *gf_array[], int num_pieces, GridFunction &merged)
2904 FiniteElementSpace *gfes = merged.FESpace();
2905 Array<int> lelem_elem, dofs;
2908 for (int i = 0; i < num_pieces; i++)
2910 FiniteElementSpace *lfes = gf_array[i]->FESpace();
2911 NURBSExtension *lext = lfes->GetMesh()->NURBSext;
2913 lext->GetElementLocalToGlobal(lelem_elem);
2915 for (int lel = 0; lel < lext->GetNE(); lel++)
2917 lfes->GetElementVDofs(lel, dofs);
2918 gf_array[i]->GetSubVector(dofs, lvec);
2920 gfes->GetElementVDofs(lelem_elem[lel], dofs);
2921 merged.SetSubVector(dofs, lvec);
2926void NURBSExtension::CheckPatches()
2928 if (Dimension() == 1 ) { return; }
2933 for (int p = 0; p < GetNP(); p++)
2935 patchTopo->GetElementEdges(p, edges, oedge);
2937 for (int i = 0; i < edges.Size(); i++)
2939 edges[i] = edge_to_knot[edges[i]];
2942 edges[i] = -1 - edges[i];
2946 if ((Dimension() == 2 &&
2947 (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
2949 (Dimension() == 3 &&
2950 (edges[0] != edges[2] || edges[0] != edges[4] ||
2951 edges[0] != edges[6] || edges[1] != edges[3] ||
2952 edges[1] != edges[5] || edges[1] != edges[7] ||
2953 edges[8] != edges[9] || edges[8] != edges[10] ||
2954 edges[8] != edges[11])))
2957 << ")\n Inconsistent edge-to-
knot mapping!\n
";
2963void NURBSExtension::CheckBdrPatches()
2968 for (int p = 0; p < GetNBP(); p++)
2970 patchTopo->GetBdrElementEdges(p, edges, oedge);
2972 for (int i = 0; i < edges.Size(); i++)
2974 edges[i] = edge_to_knot[edges[i]];
2977 edges[i] = -1 - edges[i];
2981 if ((Dimension() == 2 && (edges[0] < 0)) ||
2982 (Dimension() == 3 && (edges[0] < 0 || edges[1] < 0)))
2985 << p << ") : Bad orientation!\n
";
2991void NURBSExtension::CheckKVDirection(int p, Array <int> &kvdir)
2993 // patchTopo->GetElementEdges is not yet implemented for 1D
2994 MFEM_VERIFY(Dimension()>1, "1D not yet implemented.
");
2996 kvdir.SetSize(Dimension());
2999 Array<int> patchvert, edges, orient, edgevert;
3001 patchTopo->GetElementVertices(p, patchvert);
3003 patchTopo->GetElementEdges(p, edges, orient);
3005 // Compare the vertices of the patches with the vertices of the knotvectors of knot2dge
3006 // Based on the match the orientation will be a 1 or a -1
3007 // -1: direction is flipped
3008 // 1: direction is not flipped
3010 for (int i = 0; i < edges.Size(); i++)
3013 patchTopo->GetEdgeVertices(edges[i], edgevert);
3014 if (edgevert[0] == patchvert[0] && edgevert[1] == patchvert[1])
3019 if (edgevert[0] == patchvert[1] && edgevert[1] == patchvert[0])
3025 if (edgevert[0] == patchvert[1] && edgevert[1] == patchvert[2])
3030 if (edgevert[0] == patchvert[2] && edgevert[1] == patchvert[1])
3036 if (Dimension() == 3)
3039 for (int i = 0; i < edges.Size(); i++)
3041 patchTopo->GetEdgeVertices(edges[i], edgevert);
3043 if (edgevert[0] == patchvert[0] && edgevert[1] == patchvert[4])
3048 if (edgevert[0] == patchvert[4] && edgevert[1] == patchvert[0])
3055 MFEM_VERIFY(kvdir.Find(0) == -1, "Could not find
direction of knotvector.
");
3058void NURBSExtension::CreateComprehensiveKV()
3060 Array<int> edges, orient, kvdir;
3061 Array<int> e(Dimension());
3063 // 1D: comprehensive and unique KV are the same
3064 if (Dimension() == 1)
3066 knotVectorsCompr.SetSize(GetNKV());
3067 for (int i = 0; i < GetNKV(); i++)
3069 knotVectorsCompr[i] = new KnotVector(*(KnotVec(i)));
3073 else if (Dimension() == 2)
3075 knotVectorsCompr.SetSize(GetNP()*Dimension());
3079 else if (Dimension() == 3)
3081 knotVectorsCompr.SetSize(GetNP()*Dimension());
3087 for (int p = 0; p < GetNP(); p++)
3089 CheckKVDirection(p, kvdir);
3091 patchTopo->GetElementEdges(p, edges, orient);
3093 for (int d = 0; d < Dimension(); d++)
3095 // Indices in unique and comprehensive sets of the KnotVector
3096 int iun = edges[e[d]];
3097 int icomp = Dimension()*p+d;
3099 knotVectorsCompr[icomp] = new KnotVector(*(KnotVec(iun)));
3101 if (kvdir[d] == -1) {knotVectorsCompr[icomp]->Flip();}
3105 MFEM_VERIFY(ConsistentKVSets(), "Mismatch in KnotVectors
");
3108void NURBSExtension::UpdateUniqueKV()
3110 Array<int> e(Dimension());
3112 // 1D: comprehensive and unique KV are the same
3113 if (Dimension() == 1)
3115 for (int i = 0; i < GetNKV(); i++)
3117 *(KnotVec(i)) = *(knotVectorsCompr[i]);
3121 else if (Dimension() == 2)
3126 else if (Dimension() == 3)
3133 for (int p = 0; p < GetNP(); p++)
3135 Array<int> edges, orient, kvdir;
3137 patchTopo->GetElementEdges(p, edges, orient);
3138 CheckKVDirection(p, kvdir);
3140 for ( int d = 0; d < Dimension(); d++)
3143 if (kvdir[d] == -1) {flip = true;}
3145 // Indices in unique and comprehensive sets of the KnotVector
3146 int iun = edges[e[d]];
3147 int icomp = Dimension()*p+d;
3149 // Check if difference in order
3150 int o1 = KnotVec(iun)->GetOrder();
3151 int o2 = knotVectorsCompr[icomp]->GetOrder();
3152 int diffo = abs(o1 - o2);
3156 // Update reduced set of knotvectors
3157 *(KnotVec(iun)) = *(knotVectorsCompr[icomp]);
3159 // Give correct direction to unique knotvector.
3160 if (flip) { KnotVec(iun)->Flip(); }
3163 // Check if difference between knots
3166 if (flip) { knotVectorsCompr[icomp]->Flip(); }
3168 KnotVec(iun)->Difference(*(knotVectorsCompr[icomp]), diffknot);
3170 if (flip) { knotVectorsCompr[icomp]->Flip(); }
3172 if (diffknot.Size() > 0)
3174 // Update reduced set of knotvectors
3175 *(KnotVec(iun)) = *(knotVectorsCompr[icomp]);
3177 // Give correct direction to unique knotvector.
3178 if (flip) {KnotVec(iun)->Flip();}
3183 MFEM_VERIFY(ConsistentKVSets(), "Mismatch in KnotVectors
");
3186bool NURBSExtension::ConsistentKVSets()
3188 // patchTopo->GetElementEdges is not yet implemented for 1D
3189 MFEM_VERIFY(Dimension() > 1, "1D not yet implemented.
");
3191 Array<int> edges, orient, kvdir;
3194 Array<int> e(Dimension());
3198 if (Dimension() == 2)
3202 else if (Dimension() == 3)
3208 for (int p = 0; p < GetNP(); p++)
3210 patchTopo->GetElementEdges(p, edges, orient);
3212 CheckKVDirection(p, kvdir);
3214 for (int d = 0; d < Dimension(); d++)
3217 if (kvdir[d] == -1) {flip = true;}
3219 // Indices in unique and comprehensive sets of the KnotVector
3220 int iun = edges[e[d]];
3221 int icomp = Dimension()*p+d;
3223 // Check if KnotVectors are of equal order
3224 int o1 = KnotVec(iun)->GetOrder();
3225 int o2 = knotVectorsCompr[icomp]->GetOrder();
3226 int diffo = abs(o1 - o2);
3230 mfem::out << "\norder of knotVectorsCompr
" << d << " of patch
" << p;
3231 mfem::out << " does not agree with knotVectors
" << KnotInd(iun) << "\n
";
3235 // Check if Knotvectors have the same knots
3236 if (flip) {knotVectorsCompr[icomp]->Flip();}
3238 KnotVec(iun)->Difference(*(knotVectorsCompr[icomp]), diff);
3240 if (flip) {knotVectorsCompr[icomp]->Flip();}
3242 if (diff.Size() > 0)
3244 mfem::out << "\nknotVectorsCompr
" << d << " of patch
" << p;
3245 mfem::out << " does not agree with knotVectors
" << KnotInd(iun) << "\n
";
3253void NURBSExtension::GetPatchKnotVectors(int p, Array<KnotVector *> &kv)
3255 Array<int> edges, orient;
3257 kv.SetSize(Dimension());
3259 if (Dimension() == 1)
3261 kv[0] = knotVectorsCompr[Dimension()*p];
3263 else if (Dimension() == 2)
3265 kv[0] = knotVectorsCompr[Dimension()*p];
3266 kv[1] = knotVectorsCompr[Dimension()*p + 1];
3270 kv[0] = knotVectorsCompr[Dimension()*p];
3271 kv[1] = knotVectorsCompr[Dimension()*p + 1];
3272 kv[2] = knotVectorsCompr[Dimension()*p + 2];
3276void NURBSExtension::GetPatchKnotVectors(int p, Array<const KnotVector *> &kv)
3279 Array<int> edges, orient;
3281 kv.SetSize(Dimension());
3283 if (Dimension() == 1)
3285 kv[0] = knotVectorsCompr[Dimension()*p];
3287 else if (Dimension() == 2)
3289 kv[0] = knotVectorsCompr[Dimension()*p];
3290 kv[1] = knotVectorsCompr[Dimension()*p + 1];
3294 kv[0] = knotVectorsCompr[Dimension()*p];
3295 kv[1] = knotVectorsCompr[Dimension()*p + 1];
3296 kv[2] = knotVectorsCompr[Dimension()*p + 2];
3300void NURBSExtension::GetBdrPatchKnotVectors(int bp, Array<KnotVector *> &kv)
3305 kv.SetSize(Dimension() - 1);
3307 if (Dimension() == 2)
3309 patchTopo->GetBdrElementEdges(bp, edges, orient);
3310 kv[0] = KnotVec(edges[0]);
3312 else if (Dimension() == 3)
3314 patchTopo->GetBdrElementEdges(bp, edges, orient);
3315 kv[0] = KnotVec(edges[0]);
3316 kv[1] = KnotVec(edges[1]);
3320void NURBSExtension::GetBdrPatchKnotVectors(
3321 int bp, Array<const KnotVector *> &kv) const
3326 kv.SetSize(Dimension() - 1);
3328 if (Dimension() == 2)
3330 patchTopo->GetBdrElementEdges(bp, edges, orient);
3331 kv[0] = KnotVec(edges[0]);
3333 else if (Dimension() == 3)
3335 patchTopo->GetBdrElementEdges(bp, edges, orient);
3336 kv[0] = KnotVec(edges[0]);
3337 kv[1] = KnotVec(edges[1]);
3341void NURBSExtension::SetOrderFromOrders()
3343 MFEM_VERIFY(mOrders.Size() > 0, "");
3344 mOrder = mOrders[0];
3345 for (int i = 1; i < mOrders.Size(); i++)
3347 if (mOrders[i] != mOrder)
3349 mOrder = NURBSFECollection::VariableOrder;
3355void NURBSExtension::SetOrdersFromKnotVectors()
3357 mOrders.SetSize(NumOfKnotVectors);
3358 for (int i = 0; i < NumOfKnotVectors; i++)
3360 mOrders[i] = knotVectors[i]->GetOrder();
3362 SetOrderFromOrders();
3365void NURBSExtension::GenerateOffsets()
3367 int nv = patchTopo->GetNV();
3368 int ne = patchTopo->GetNEdges();
3369 int nf = patchTopo->GetNFaces();
3370 int np = patchTopo->GetNE();
3371 int meshCounter, spaceCounter, dim = Dimension();
3376 v_meshOffsets.SetSize(nv);
3377 e_meshOffsets.SetSize(ne);
3378 f_meshOffsets.SetSize(nf);
3379 p_meshOffsets.SetSize(np);
3381 v_spaceOffsets.SetSize(nv);
3382 e_spaceOffsets.SetSize(ne);
3383 f_spaceOffsets.SetSize(nf);
3384 p_spaceOffsets.SetSize(np);
3386 // Get vertex offsets
3387 for (meshCounter = 0; meshCounter < nv; meshCounter++)
3389 v_meshOffsets[meshCounter] = meshCounter;
3390 v_spaceOffsets[meshCounter] = meshCounter;
3392 spaceCounter = meshCounter;
3395 for (int e = 0; e < ne; e++)
3397 e_meshOffsets[e] = meshCounter;
3398 e_spaceOffsets[e] = spaceCounter;
3399 meshCounter += KnotVec(e)->GetNE() - 1;
3400 spaceCounter += KnotVec(e)->GetNCP() - 2;
3404 for (int f = 0; f < nf; f++)
3406 f_meshOffsets[f] = meshCounter;
3407 f_spaceOffsets[f] = spaceCounter;
3409 patchTopo->GetFaceEdges(f, edges, orient);
3412 (KnotVec(edges[0])->GetNE() - 1) *
3413 (KnotVec(edges[1])->GetNE() - 1);
3415 (KnotVec(edges[0])->GetNCP() - 2) *
3416 (KnotVec(edges[1])->GetNCP() - 2);
3419 // Get patch offsets
3420 for (int p = 0; p < np; p++)
3422 p_meshOffsets[p] = meshCounter;
3423 p_spaceOffsets[p] = spaceCounter;
3427 meshCounter += KnotVec(0)->GetNE() - 1;
3428 spaceCounter += KnotVec(0)->GetNCP() - 2;
3432 patchTopo->GetElementEdges(p, edges, orient);
3434 (KnotVec(edges[0])->GetNE() - 1) *
3435 (KnotVec(edges[1])->GetNE() - 1);
3437 (KnotVec(edges[0])->GetNCP() - 2) *
3438 (KnotVec(edges[1])->GetNCP() - 2);
3442 patchTopo->GetElementEdges(p, edges, orient);
3444 (KnotVec(edges[0])->GetNE() - 1) *
3445 (KnotVec(edges[3])->GetNE() - 1) *
3446 (KnotVec(edges[8])->GetNE() - 1);
3448 (KnotVec(edges[0])->GetNCP() - 2) *
3449 (KnotVec(edges[3])->GetNCP() - 2) *
3450 (KnotVec(edges[8])->GetNCP() - 2);
3453 NumOfVertices = meshCounter;
3454 NumOfDofs = spaceCounter;
3457void NURBSExtension::CountElements()
3459 int dim = Dimension();
3460 Array<const KnotVector *> kv(dim);
3463 for (int p = 0; p < GetNP(); p++)
3465 GetPatchKnotVectors(p, kv);
3467 int ne = kv[0]->GetNE();
3468 for (int d = 1; d < dim; d++)
3470 ne *= kv[d]->GetNE();
3473 NumOfElements += ne;
3477void NURBSExtension::CountBdrElements()
3479 int dim = Dimension() - 1;
3480 Array<KnotVector *> kv(dim);
3482 NumOfBdrElements = 0;
3483 for (int p = 0; p < GetNBP(); p++)
3485 GetBdrPatchKnotVectors(p, kv);
3488 for (int d = 0; d < dim; d++)
3490 ne *= kv[d]->GetNE();
3493 NumOfBdrElements += ne;
3497void NURBSExtension::GetElementTopo(Array<Element *> &elements) const
3499 elements.SetSize(GetNE());
3501 if (Dimension() == 1)
3503 Get1DElementTopo(elements);
3505 else if (Dimension() == 2)
3507 Get2DElementTopo(elements);
3511 Get3DElementTopo(elements);
3515void NURBSExtension::Get1DElementTopo(Array<Element *> &elements) const
3520 NURBSPatchMap p2g(this);
3521 const KnotVector *kv[1];
3523 for (int p = 0; p < GetNP(); p++)
3525 p2g.SetPatchVertexMap(p, kv);
3528 int patch_attr = patchTopo->GetAttribute(p);
3530 for (int i = 0; i < nx; i++)
3534 ind[0] = activeVert[p2g(i)];
3535 ind[1] = activeVert[p2g(i+1)];
3537 elements[el] = new Segment(ind, patch_attr);
3545void NURBSExtension::Get2DElementTopo(Array<Element *> &elements) const
3550 NURBSPatchMap p2g(this);
3551 const KnotVector *kv[2];
3553 for (int p = 0; p < GetNP(); p++)
3555 p2g.SetPatchVertexMap(p, kv);
3559 int patch_attr = patchTopo->GetAttribute(p);
3561 for (int j = 0; j < ny; j++)
3563 for (int i = 0; i < nx; i++)
3567 ind[0] = activeVert[p2g(i, j )];
3568 ind[1] = activeVert[p2g(i+1,j )];
3569 ind[2] = activeVert[p2g(i+1,j+1)];
3570 ind[3] = activeVert[p2g(i, j+1)];
3572 elements[el] = new Quadrilateral(ind, patch_attr);
3581void NURBSExtension::Get3DElementTopo(Array<Element *> &elements) const
3586 NURBSPatchMap p2g(this);
3587 const KnotVector *kv[3];
3589 for (int p = 0; p < GetNP(); p++)
3591 p2g.SetPatchVertexMap(p, kv);
3596 int patch_attr = patchTopo->GetAttribute(p);
3598 for (int k = 0; k < nz; k++)
3600 for (int j = 0; j < ny; j++)
3602 for (int i = 0; i < nx; i++)
3606 ind[0] = activeVert[p2g(i, j, k)];
3607 ind[1] = activeVert[p2g(i+1,j, k)];
3608 ind[2] = activeVert[p2g(i+1,j+1,k)];
3609 ind[3] = activeVert[p2g(i, j+1,k)];
3611 ind[4] = activeVert[p2g(i, j, k+1)];
3612 ind[5] = activeVert[p2g(i+1,j, k+1)];
3613 ind[6] = activeVert[p2g(i+1,j+1,k+1)];
3614 ind[7] = activeVert[p2g(i, j+1,k+1)];
3616 elements[el] = new Hexahedron(ind, patch_attr);
3626void NURBSExtension::GetBdrElementTopo(Array<Element *> &boundary) const
3628 boundary.SetSize(GetNBE());
3630 if (Dimension() == 1)
3632 Get1DBdrElementTopo(boundary);
3634 else if (Dimension() == 2)
3636 Get2DBdrElementTopo(boundary);
3640 Get3DBdrElementTopo(boundary);
3644void NURBSExtension::Get1DBdrElementTopo(Array<Element *> &boundary) const
3648 NURBSPatchMap p2g(this);
3649 const KnotVector *kv[1];
3652 for (int b = 0; b < GetNBP(); b++)
3654 p2g.SetBdrPatchVertexMap(b, kv, okv);
3655 int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
3657 if (activeBdrElem[g_be])
3659 ind[0] = activeVert[p2g[0]];
3660 boundary[l_be] = new Point(ind, bdr_patch_attr);
3667void NURBSExtension::Get2DBdrElementTopo(Array<Element *> &boundary) const
3671 NURBSPatchMap p2g(this);
3672 const KnotVector *kv[1];
3675 for (int b = 0; b < GetNBP(); b++)
3677 p2g.SetBdrPatchVertexMap(b, kv, okv);
3680 int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
3682 for (int i = 0; i < nx; i++)
3684 if (activeBdrElem[g_be])
3686 int i_ = (okv[0] >= 0) ? i : (nx - 1 - i);
3687 ind[0] = activeVert[p2g[i_ ]];
3688 ind[1] = activeVert[p2g[i_+1]];
3690 boundary[l_be] = new Segment(ind, bdr_patch_attr);
3698void NURBSExtension::Get3DBdrElementTopo(Array<Element *> &boundary) const
3702 NURBSPatchMap p2g(this);
3703 const KnotVector *kv[2];
3706 for (int b = 0; b < GetNBP(); b++)
3708 p2g.SetBdrPatchVertexMap(b, kv, okv);
3712 int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
3714 for (int j = 0; j < ny; j++)
3716 int j_ = (okv[1] >= 0) ? j : (ny - 1 - j);
3717 for (int i = 0; i < nx; i++)
3719 if (activeBdrElem[g_be])
3721 int i_ = (okv[0] >= 0) ? i : (nx - 1 - i);
3722 ind[0] = activeVert[p2g(i_, j_ )];
3723 ind[1] = activeVert[p2g(i_+1,j_ )];
3724 ind[2] = activeVert[p2g(i_+1,j_+1)];
3725 ind[3] = activeVert[p2g(i_, j_+1)];
3727 boundary[l_be] = new Quadrilateral(ind, bdr_patch_attr);
3736void NURBSExtension::GenerateElementDofTable()
3738 activeDof.SetSize(GetNTotalDof());
3741 if (Dimension() == 1)
3743 Generate1DElementDofTable();
3745 else if (Dimension() == 2)
3747 Generate2DElementDofTable();
3751 Generate3DElementDofTable();
3754 SetPatchToElements();
3756 NumOfActiveDofs = 0;
3757 for (int d = 0; d < GetNTotalDof(); d++)
3761 activeDof[d] = NumOfActiveDofs;
3764 int *dof = el_dof->GetJ();
3765 int ndof = el_dof->Size_of_connections();
3766 for (int i = 0; i < ndof; i++)
3768 dof[i] = activeDof[dof[i]] - 1;
3772void NURBSExtension::Generate1DElementDofTable()
3776 const KnotVector *kv[2];
3777 NURBSPatchMap p2g(this);
3779 Array<Connection> el_dof_list;
3780 el_to_patch.SetSize(NumOfActiveElems);
3781 el_to_IJK.SetSize(NumOfActiveElems, 2);
3783 for (int p = 0; p < GetNP(); p++)
3785 p2g.SetPatchDofMap(p, kv);
3788 const int ord0 = kv[0]->GetOrder();
3789 for (int i = 0; i < kv[0]->GetNKS(); i++)
3791 if (kv[0]->isElement(i))
3795 Connection conn(el,0);
3796 for (int ii = 0; ii <= ord0; ii++)
3798 conn.to = DofMap(p2g(i+ii));
3799 activeDof[conn.to] = 1;
3800 el_dof_list.Append(conn);
3802 el_to_patch[el] = p;
3803 el_to_IJK(el,0) = i;
3811 // We must NOT sort el_dof_list in this case.
3812 el_dof = new Table(NumOfActiveElems, el_dof_list);
3815void NURBSExtension::Generate2DElementDofTable()
3819 const KnotVector *kv[2];
3820 NURBSPatchMap p2g(this);
3822 Array<Connection> el_dof_list;
3823 el_to_patch.SetSize(NumOfActiveElems);
3824 el_to_IJK.SetSize(NumOfActiveElems, 2);
3826 for (int p = 0; p < GetNP(); p++)
3828 p2g.SetPatchDofMap(p, kv);
3831 const int ord0 = kv[0]->GetOrder();
3832 const int ord1 = kv[1]->GetOrder();
3833 for (int j = 0; j < kv[1]->GetNKS(); j++)
3835 if (kv[1]->isElement(j))
3837 for (int i = 0; i < kv[0]->GetNKS(); i++)
3839 if (kv[0]->isElement(i))
3843 Connection conn(el,0);
3844 for (int jj = 0; jj <= ord1; jj++)
3846 for (int ii = 0; ii <= ord0; ii++)
3848 conn.to = DofMap(p2g(i+ii,j+jj));
3849 activeDof[conn.to] = 1;
3850 el_dof_list.Append(conn);
3853 el_to_patch[el] = p;
3854 el_to_IJK(el,0) = i;
3855 el_to_IJK(el,1) = j;
3865 // We must NOT sort el_dof_list in this case.
3866 el_dof = new Table(NumOfActiveElems, el_dof_list);
3869void NURBSExtension::Generate3DElementDofTable()
3873 const KnotVector *kv[3];
3874 NURBSPatchMap p2g(this);
3876 Array<Connection> el_dof_list;
3877 el_to_patch.SetSize(NumOfActiveElems);
3878 el_to_IJK.SetSize(NumOfActiveElems, 3);
3880 for (int p = 0; p < GetNP(); p++)
3882 p2g.SetPatchDofMap(p, kv);
3885 const int ord0 = kv[0]->GetOrder();
3886 const int ord1 = kv[1]->GetOrder();
3887 const int ord2 = kv[2]->GetOrder();
3888 for (int k = 0; k < kv[2]->GetNKS(); k++)
3890 if (kv[2]->isElement(k))
3892 for (int j = 0; j < kv[1]->GetNKS(); j++)
3894 if (kv[1]->isElement(j))
3896 for (int i = 0; i < kv[0]->GetNKS(); i++)
3898 if (kv[0]->isElement(i))
3902 Connection conn(el,0);
3903 for (int kk = 0; kk <= ord2; kk++)
3905 for (int jj = 0; jj <= ord1; jj++)
3907 for (int ii = 0; ii <= ord0; ii++)
3909 conn.to = DofMap(p2g(i+ii, j+jj, k+kk));
3910 activeDof[conn.to] = 1;
3911 el_dof_list.Append(conn);
3916 el_to_patch[el] = p;
3917 el_to_IJK(el,0) = i;
3918 el_to_IJK(el,1) = j;
3919 el_to_IJK(el,2) = k;
3931 // We must NOT sort el_dof_list in this case.
3932 el_dof = new Table(NumOfActiveElems, el_dof_list);
3935void NURBSExtension::GetPatchDofs(const int patch, Array<int> &dofs)
3937 const KnotVector *kv[3];
3938 NURBSPatchMap p2g(this);
3940 p2g.SetPatchDofMap(patch, kv);
3942 if (Dimension() == 1)
3944 const int nx = kv[0]->GetNCP();
3947 for (int i=0; i<nx; ++i)
3949 dofs[i] = DofMap(p2g(i));
3952 else if (Dimension() == 2)
3954 const int nx = kv[0]->GetNCP();
3955 const int ny = kv[1]->GetNCP();
3956 dofs.SetSize(nx * ny);
3958 for (int j=0; j<ny; ++j)
3959 for (int i=0; i<nx; ++i)
3961 dofs[i + (nx * j)] = DofMap(p2g(i, j));
3964 else if (Dimension() == 3)
3966 const int nx = kv[0]->GetNCP();
3967 const int ny = kv[1]->GetNCP();
3968 const int nz = kv[2]->GetNCP();
3969 dofs.SetSize(nx * ny * nz);
3971 for (int k=0; k<nz; ++k)
3972 for (int j=0; j<ny; ++j)
3973 for (int i=0; i<nx; ++i)
3975 dofs[i + (nx * (j + (k * ny)))] = DofMap(p2g(i, j, k));
3980 MFEM_ABORT("Only 1D/2D/3D supported currently in
NURBSExtension::GetPatchDofs
");
3984void NURBSExtension::GenerateBdrElementDofTable()
3986 if (Dimension() == 1)
3988 Generate1DBdrElementDofTable();
3990 else if (Dimension() == 2)
3992 Generate2DBdrElementDofTable();
3996 Generate3DBdrElementDofTable();
3999 SetPatchToBdrElements();
4001 int *dof = bel_dof->GetJ();
4002 int ndof = bel_dof->Size_of_connections();
4003 for (int i = 0; i < ndof; i++)
4008 dof[i] = -1 - (activeDof[-1-idx] - 1);
4009 dof[i] = -activeDof[-1-idx];
4013 dof[i] = activeDof[idx] - 1;
4018void NURBSExtension::Generate1DBdrElementDofTable()
4021 int lbe = 0, okv[1];
4022 const KnotVector *kv[1];
4023 NURBSPatchMap p2g(this);
4025 Array<Connection> bel_dof_list;
4026 bel_to_patch.SetSize(NumOfActiveBdrElems);
4027 bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
4029 for (int b = 0; b < GetNBP(); b++)
4031 p2g.SetBdrPatchDofMap(b, kv, okv);
4033 if (activeBdrElem[gbe])
4035 Connection conn(lbe,0);
4036 conn.to = DofMap(p2g[0]);
4037 bel_dof_list.Append(conn);
4038 bel_to_patch[lbe] = b;
4039 bel_to_IJK(lbe,0) = 0;
4044 // We must NOT sort bel_dof_list in this case.
4045 bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
4048void NURBSExtension::Generate2DBdrElementDofTable()
4051 int lbe = 0, okv[1];
4052 const KnotVector *kv[1];
4053 NURBSPatchMap p2g(this);
4055 Array<Connection> bel_dof_list;
4056 bel_to_patch.SetSize(NumOfActiveBdrElems);
4057 bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
4059 for (int b = 0; b < GetNBP(); b++)
4061 p2g.SetBdrPatchDofMap(b, kv, okv);
4062 const int nx = p2g.nx(); // NCP-1
4064 const int nks0 = kv[0]->GetNKS();
4065 const int ord0 = kv[0]->GetOrder();
4067 bool add_dofs = true;
4070 if (mode == Mode::H_DIV)
4072 int fn = patchTopo->GetBdrElementFaceIndex(b);
4073 if (ord0 == mOrders.Max()) { add_dofs = false; }
4074 if (fn == 0) { s = -1; }
4075 if (fn == 2) { s = -1; }
4077 else if (mode == Mode::H_CURL)
4079 if (ord0 == mOrders.Max()) { add_dofs = false; }
4082 for (int i = 0; i < nks0; i++)
4084 if (kv[0]->isElement(i))
4086 if (activeBdrElem[gbe])
4088 Connection conn(lbe,0);
4091 for (int ii = 0; ii <= ord0; ii++)
4093 conn.to = DofMap(p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)]);
4094 if (s == -1) { conn.to = -1 -conn.to; }
4095 bel_dof_list.Append(conn);
4098 bel_to_patch[lbe] = b;
4099 bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
4106 // We must NOT sort bel_dof_list in this case.
4107 bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
4111void NURBSExtension::Generate3DBdrElementDofTable()
4114 int lbe = 0, okv[2];
4115 const KnotVector *kv[2];
4116 NURBSPatchMap p2g(this);
4118 Array<Connection> bel_dof_list;
4119 bel_to_patch.SetSize(NumOfActiveBdrElems);
4120 bel_to_IJK.SetSize(NumOfActiveBdrElems, 2);
4122 for (int b = 0; b < GetNBP(); b++)
4124 p2g.SetBdrPatchDofMap(b, kv, okv);
4125 const int nx = p2g.nx(); // NCP0-1
4126 const int ny = p2g.ny(); // NCP1-1
4129 const int nks0 = kv[0]->GetNKS();
4130 const int ord0 = kv[0]->GetOrder();
4131 const int nks1 = kv[1]->GetNKS();
4132 const int ord1 = kv[1]->GetOrder();
4134 // Check if dofs are actually defined on boundary
4135 bool add_dofs = true;
4138 if (mode == Mode::H_DIV)
4140 int fn = patchTopo->GetBdrElementFaceIndex(b);
4141 if (ord0 != ord1) { add_dofs = false; }
4142 if (fn == 4) { s = -1; }
4143 if (fn == 1) { s = -1; }
4144 if (fn == 0) { s = -1; }
4146 else if (mode == Mode::H_CURL)
4148 if (ord0 == ord1) { add_dofs = false; }
4152 for (int j = 0; j < nks1; j++)
4154 if (kv[1]->isElement(j))
4156 for (int i = 0; i < nks0; i++)
4158 if (kv[0]->isElement(i))
4160 if (activeBdrElem[gbe])
4162 Connection conn(lbe,0);
4165 for (int jj = 0; jj <= ord1; jj++)
4167 const int jj_ = (okv[1] >= 0) ? (j+jj) : (ny-j-jj);
4168 for (int ii = 0; ii <= ord0; ii++)
4170 const int ii_ = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
4171 conn.to = DofMap(p2g(ii_, jj_));
4172 if (s == -1) { conn.to = -1 -conn.to; }
4173 bel_dof_list.Append(conn);
4177 bel_to_patch[lbe] = b;
4178 bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
4179 bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
4188 // We must NOT sort bel_dof_list in this case.
4189 bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
4192void NURBSExtension::GetVertexLocalToGlobal(Array<int> &lvert_vert)
4194 lvert_vert.SetSize(GetNV());
4195 for (int gv = 0; gv < GetGNV(); gv++)
4196 if (activeVert[gv] >= 0)
4198 lvert_vert[activeVert[gv]] = gv;
4202void NURBSExtension::GetElementLocalToGlobal(Array<int> &lelem_elem)
4204 lelem_elem.SetSize(GetNE());
4205 for (int le = 0, ge = 0; ge < GetGNE(); ge++)
4208 lelem_elem[le++] = ge;
4212void NURBSExtension::LoadFE(int i, const FiniteElement *FE) const
4214 const NURBSFiniteElement *NURBSFE =
4215 dynamic_cast<const NURBSFiniteElement *>(FE);
4217 if (NURBSFE->GetElement() != i)
4220 NURBSFE->SetIJK(el_to_IJK.GetRow(i));
4221 if (el_to_patch[i] != NURBSFE->GetPatch())
4223 GetPatchKnotVectors(el_to_patch[i], NURBSFE->KnotVectors());
4224 NURBSFE->SetPatch(el_to_patch[i]);
4225 NURBSFE->SetOrder();
4227 el_dof->GetRow(i, dofs);
4228 weights.GetSubVector(dofs, NURBSFE->Weights());
4229 NURBSFE->SetElement(i);
4233void NURBSExtension::LoadBE(int i, const FiniteElement *BE) const
4235 if (Dimension() == 1) { return; }
4237 const NURBSFiniteElement *NURBSFE =
4238 dynamic_cast<const NURBSFiniteElement *>(BE);
4240 if (NURBSFE->GetElement() != i)
4243 NURBSFE->SetIJK(bel_to_IJK.GetRow(i));
4244 if (bel_to_patch[i] != NURBSFE->GetPatch())
4246 GetBdrPatchKnotVectors(bel_to_patch[i], NURBSFE->KnotVectors());
4247 NURBSFE->SetPatch(bel_to_patch[i]);
4248 NURBSFE->SetOrder();
4250 bel_dof->GetRow(i, dofs);
4251 weights.GetSubVector(dofs, NURBSFE->Weights());
4252 NURBSFE->SetElement(i);
4256void NURBSExtension::ConvertToPatches(const Vector &Nodes)
4261 if (patches.Size() == 0)
4263 GetPatchNets(Nodes, Dimension());
4267void NURBSExtension::SetCoordsFromPatches(Vector &Nodes)
4269 if (patches.Size() == 0) { return; }
4271 SetSolutionVector(Nodes, Dimension());
4275void NURBSExtension::SetKnotsFromPatches()
4277 if (patches.Size() == 0)
4280 " No patches available!
");
4283 Array<KnotVector *> kv;
4285 for (int p = 0; p < patches.Size(); p++)
4287 GetPatchKnotVectors(p, kv);
4289 for (int i = 0; i < kv.Size(); i++)
4291 *kv[i] = *patches[p]->GetKV(i);
4296 SetOrdersFromKnotVectors();
4302 // all elements must be active
4303 NumOfActiveElems = NumOfElements;
4304 activeElem.SetSize(NumOfElements);
4307 GenerateActiveVertices();
4309 GenerateElementDofTable();
4310 GenerateActiveBdrElems();
4311 GenerateBdrElementDofTable();
4313 ConnectBoundaries();
4316void NURBSExtension::LoadSolution(std::istream &input, GridFunction &sol) const
4318 const FiniteElementSpace *fes = sol.FESpace();
4319 MFEM_VERIFY(fes->GetNURBSext() == this, "");
4321 sol.SetSize(fes->GetVSize());
4323 Array<const KnotVector *> kv(Dimension());
4324 NURBSPatchMap p2g(this);
4325 const int vdim = fes->GetVDim();
4327 for (int p = 0; p < GetNP(); p++)
4329 skip_comment_lines(input, '#');
4331 p2g.SetPatchDofMap(p, kv);
4332 const int nx = kv[0]->GetNCP();
4333 const int ny = kv[1]->GetNCP();
4334 const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
4335 for (int k = 0; k < nz; k++)
4337 for (int j = 0; j < ny; j++)
4339 for (int i = 0; i < nx; i++)
4341 const int ll = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
4342 const int l = DofMap(ll);
4343 for (int vd = 0; vd < vdim; vd++)
4345 input >> sol(fes->DofToVDof(l,vd));
4353void NURBSExtension::PrintSolution(const GridFunction &sol, std::ostream &os)
4356 const FiniteElementSpace *fes = sol.FESpace();
4357 MFEM_VERIFY(fes->GetNURBSext() == this, "");
4359 Array<const KnotVector *> kv(Dimension());
4360 NURBSPatchMap p2g(this);
4361 const int vdim = fes->GetVDim();
4363 for (int p = 0; p < GetNP(); p++)
4365 os << "\n# patch
" << p << "\n\n
";
4367 p2g.SetPatchDofMap(p, kv);
4368 const int nx = kv[0]->GetNCP();
4369 const int ny = kv[1]->GetNCP();
4370 const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
4371 for (int k = 0; k < nz; k++)
4373 for (int j = 0; j < ny; j++)
4375 for (int i = 0; i < nx; i++)
4377 const int ll = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
4378 const int l = DofMap(ll);
4379 os << sol(fes->DofToVDof(l,0));
4380 for (int vd = 1; vd < vdim; vd++)
4382 os << ' ' << sol(fes->DofToVDof(l,vd));
4391void NURBSExtension::DegreeElevate(int rel_degree, int degree)
4393 for (int p = 0; p < patches.Size(); p++)
4395 for (int dir = 0; dir < patches[p]->GetNKV(); dir++)
4397 int oldd = patches[p]->GetKV(dir)->GetOrder();
4398 int newd = std::min(oldd + rel_degree, degree);
4401 patches[p]->DegreeElevate(dir, newd - oldd);
4407NURBSExtension* NURBSExtension::GetDivExtension(int component)
4413 "only works for single patch NURBS meshes
");
4416 Array<int> newOrders = GetOrders();
4417 newOrders[component] += 1;
4419 return new NURBSExtension(this, newOrders, Mode::H_DIV);
4422NURBSExtension* NURBSExtension::GetCurlExtension(int component)
4428 "only works for single patch NURBS meshes
");
4431 Array<int> newOrders = GetOrders();
4432 for (int c = 0; c < newOrders.Size(); c++) { newOrders[c]++; }
4433 newOrders[component] -= 1;
4435 return new NURBSExtension(this, newOrders, Mode::H_CURL);
4439void NURBSExtension::UniformRefinement(Array<int> const& rf)
4441 for (int p = 0; p < patches.Size(); p++)
4443 patches[p]->UniformRefinement(rf);
4447void NURBSExtension::UniformRefinement(int rf)
4449 Array<int> rf_array(Dimension());
4451 UniformRefinement(rf_array);
4454void NURBSExtension::Coarsen(Array<int> const& cf, real_t tol)
4456 // First, mark all knot vectors on all patches as not coarse. This prevents
4457 // coarsening the same knot vector twice.
4458 for (int p = 0; p < patches.Size(); p++)
4460 patches[p]->SetKnotVectorsCoarse(false);
4463 for (int p = 0; p < patches.Size(); p++)
4465 patches[p]->Coarsen(cf, tol);
4469void NURBSExtension::Coarsen(int cf, real_t tol)
4471 Array<int> cf_array(Dimension());
4473 Coarsen(cf_array, tol);
4476void NURBSExtension::GetCoarseningFactors(Array<int> & f) const
4479 for (auto patch : patches)
4482 patch->GetCoarseningFactors(pf);
4485 f = pf; // Initialize
4489 MFEM_VERIFY(f.Size() == pf.Size(), "");
4490 for (int i=0; i<f.Size(); ++i)
4492 MFEM_VERIFY(f[i] == pf[i] || f[i] == 1 || pf[i] == 1,
4493 "Inconsistent patch coarsening factors
");
4494 if (f[i] == 1 && pf[i] != 1)
4503void NURBSExtension::KnotInsert(Array<KnotVector *> &kv)
4509 Array<KnotVector *> pkv(Dimension());
4511 for (int p = 0; p < patches.Size(); p++)
4515 pkv[0] = kv[KnotInd(p)];
4517 else if (Dimension()==2)
4519 patchTopo->GetElementEdges(p, edges, orient);
4520 pkv[0] = kv[KnotInd(edges[0])];
4521 pkv[1] = kv[KnotInd(edges[1])];
4523 else if (Dimension()==3)
4525 patchTopo->GetElementEdges(p, edges, orient);
4526 pkv[0] = kv[KnotInd(edges[0])];
4527 pkv[1] = kv[KnotInd(edges[3])];
4528 pkv[2] = kv[KnotInd(edges[8])];
4531 // Check whether inserted knots should be flipped before inserting.
4532 // Knotvectors are stored in a different array pkvc such that the original
4533 // knots which are inserted are not changed.
4534 // We need those knots for multiple patches so they have to remain original
4535 CheckKVDirection(p, kvdir);
4537 Array<KnotVector *> pkvc(Dimension());
4538 for (int d = 0; d < Dimension(); d++)
4540 pkvc[d] = new KnotVector(*(pkv[d]));
4548 patches[p]->KnotInsert(pkvc);
4549 for (int d = 0; d < Dimension(); d++) { delete pkvc[d]; }
4553void NURBSExtension::KnotInsert(Array<Vector *> &kv)
4559 Array<Vector *> pkv(Dimension());
4561 for (int p = 0; p < patches.Size(); p++)
4565 pkv[0] = kv[KnotInd(p)];
4567 else if (Dimension()==2)
4569 patchTopo->GetElementEdges(p, edges, orient);
4570 pkv[0] = kv[KnotInd(edges[0])];
4571 pkv[1] = kv[KnotInd(edges[1])];
4573 else if (Dimension()==3)
4575 patchTopo->GetElementEdges(p, edges, orient);
4576 pkv[0] = kv[KnotInd(edges[0])];
4577 pkv[1] = kv[KnotInd(edges[3])];
4578 pkv[2] = kv[KnotInd(edges[8])];
4581 // Check whether inserted knots should be flipped before inserting.
4582 // Knotvectors are stored in a different array pkvc such that the original
4583 // knots which are inserted are not changed.
4584 CheckKVDirection(p, kvdir);
4586 Array<Vector *> pkvc(Dimension());
4587 for (int d = 0; d < Dimension(); d++)
4589 pkvc[d] = new Vector(*(pkv[d]));
4593 // Find flip point, for knotvectors that do not have the domain [0:1]
4594 KnotVector *kva = knotVectorsCompr[Dimension()*p+d];
4595 real_t apb = (*kva)[0] + (*kva)[kva->Size()-1];
4598 int size = pkvc[d]->Size();
4599 int ns = ceil(size/2.0);
4600 for (int j = 0; j < ns; j++)
4602 real_t tmp = apb - pkvc[d]->Elem(j);
4603 pkvc[d]->Elem(j) = apb - pkvc[d]->Elem(size-1-j);
4604 pkvc[d]->Elem(size-1-j) = tmp;
4609 patches[p]->KnotInsert(pkvc);
4611 for (int i = 0; i < Dimension(); i++) { delete pkvc[i]; }
4615void NURBSExtension::KnotRemove(Array<Vector *> &kv, real_t tol)
4621 Array<Vector *> pkv(Dimension());
4623 for (int p = 0; p < patches.Size(); p++)
4627 pkv[0] = kv[KnotInd(p)];
4629 else if (Dimension()==2)
4631 patchTopo->GetElementEdges(p, edges, orient);
4632 pkv[0] = kv[KnotInd(edges[0])];
4633 pkv[1] = kv[KnotInd(edges[1])];
4635 else if (Dimension()==3)
4637 patchTopo->GetElementEdges(p, edges, orient);
4638 pkv[0] = kv[KnotInd(edges[0])];
4639 pkv[1] = kv[KnotInd(edges[3])];
4640 pkv[2] = kv[KnotInd(edges[8])];
4643 // Check whether knots should be flipped before removing.
4644 CheckKVDirection(p, kvdir);
4646 Array<Vector *> pkvc(Dimension());
4647 for (int d = 0; d < Dimension(); d++)
4649 pkvc[d] = new Vector(*(pkv[d]));
4653 // Find flip point, for knotvectors that do not have the domain [0:1]
4654 KnotVector *kva = knotVectorsCompr[Dimension()*p+d];
4655 real_t apb = (*kva)[0] + (*kva)[kva->Size()-1];
4658 int size = pkvc[d]->Size();
4659 int ns = ceil(size/2.0);
4660 for (int j = 0; j < ns; j++)
4662 real_t tmp = apb - pkvc[d]->Elem(j);
4663 pkvc[d]->Elem(j) = apb - pkvc[d]->Elem(size-1-j);
4664 pkvc[d]->Elem(size-1-j) = tmp;
4669 patches[p]->KnotRemove(pkvc, tol);
4671 for (int i = 0; i < Dimension(); i++) { delete pkvc[i]; }
4675void NURBSExtension::GetPatchNets(const Vector &coords, int vdim)
4677 if (Dimension() == 1)
4679 Get1DPatchNets(coords, vdim);
4681 else if (Dimension() == 2)
4683 Get2DPatchNets(coords, vdim);
4687 Get3DPatchNets(coords, vdim);
4691void NURBSExtension::Get1DPatchNets(const Vector &coords, int vdim)
4693 Array<const KnotVector *> kv(1);
4694 NURBSPatchMap p2g(this);
4696 patches.SetSize(GetNP());
4697 for (int p = 0; p < GetNP(); p++)
4699 p2g.SetPatchDofMap(p, kv);
4700 patches[p] = new NURBSPatch(kv, vdim+1);
4701 NURBSPatch &Patch = *patches[p];
4703 for (int i = 0; i < kv[0]->GetNCP(); i++)
4705 const int l = DofMap(p2g(i));
4706 for (int d = 0; d < vdim; d++)
4708 Patch(i,d) = coords(l*vdim + d)*weights(l);
4710 Patch(i,vdim) = weights(l);
4716void NURBSExtension::Get2DPatchNets(const Vector &coords, int vdim)
4718 Array<const KnotVector *> kv(2);
4719 NURBSPatchMap p2g(this);
4721 patches.SetSize(GetNP());
4722 for (int p = 0; p < GetNP(); p++)
4724 p2g.SetPatchDofMap(p, kv);
4725 patches[p] = new NURBSPatch(kv, vdim+1);
4726 NURBSPatch &Patch = *patches[p];
4728 for (int j = 0; j < kv[1]->GetNCP(); j++)
4730 for (int i = 0; i < kv[0]->GetNCP(); i++)
4732 const int l = DofMap(p2g(i,j));
4733 for (int d = 0; d < vdim; d++)
4735 Patch(i,j,d) = coords(l*vdim + d)*weights(l);
4737 Patch(i,j,vdim) = weights(l);
4743void NURBSExtension::Get3DPatchNets(const Vector &coords, int vdim)
4745 Array<const KnotVector *> kv(3);
4746 NURBSPatchMap p2g(this);
4748 patches.SetSize(GetNP());
4749 for (int p = 0; p < GetNP(); p++)
4751 p2g.SetPatchDofMap(p, kv);
4752 patches[p] = new NURBSPatch(kv, vdim+1);
4753 NURBSPatch &Patch = *patches[p];
4755 for (int k = 0; k < kv[2]->GetNCP(); k++)
4757 for (int j = 0; j < kv[1]->GetNCP(); j++)
4759 for (int i = 0; i < kv[0]->GetNCP(); i++)
4761 const int l = DofMap(p2g(i,j,k));
4762 for (int d = 0; d < vdim; d++)
4764 Patch(i,j,k,d) = coords(l*vdim + d)*weights(l);
4766 Patch(i,j,k,vdim) = weights(l);
4773void NURBSExtension::SetSolutionVector(Vector &coords, int vdim)
4775 if (Dimension() == 1)
4777 Set1DSolutionVector(coords, vdim);
4779 else if (Dimension() == 2)
4781 Set2DSolutionVector(coords, vdim);
4785 Set3DSolutionVector(coords, vdim);
4789void NURBSExtension::Set1DSolutionVector(Vector &coords, int vdim)
4791 Array<const KnotVector *> kv(1);
4792 NURBSPatchMap p2g(this);
4794 weights.SetSize(GetNDof());
4795 for (int p = 0; p < GetNP(); p++)
4797 p2g.SetPatchDofMap(p, kv);
4798 NURBSPatch &patch = *patches[p];
4799 MFEM_ASSERT(vdim+1 == patch.GetNC(), "");
4801 for (int i = 0; i < kv[0]->GetNCP(); i++)
4803 const int l = p2g(i);
4804 for (int d = 0; d < vdim; d++)
4806 coords(l*vdim + d) = patch(i,d)/patch(i,vdim);
4808 weights(l) = patch(i,vdim);
4816void NURBSExtension::Set2DSolutionVector(Vector &coords, int vdim)
4818 Array<const KnotVector *> kv(2);
4819 NURBSPatchMap p2g(this);
4821 weights.SetSize(GetNDof());
4822 for (int p = 0; p < GetNP(); p++)
4824 p2g.SetPatchDofMap(p, kv);
4825 NURBSPatch &patch = *patches[p];
4826 MFEM_ASSERT(vdim+1 == patch.GetNC(), "");
4828 for (int j = 0; j < kv[1]->GetNCP(); j++)
4830 for (int i = 0; i < kv[0]->GetNCP(); i++)
4832 const int l = p2g(i,j);
4833 for (int d = 0; d < vdim; d++)
4835 coords(l*vdim + d) = patch(i,j,d)/patch(i,j,vdim);
4837 weights(l) = patch(i,j,vdim);
4844void NURBSExtension::Set3DSolutionVector(Vector &coords, int vdim)
4846 Array<const KnotVector *> kv(3);
4847 NURBSPatchMap p2g(this);
4849 weights.SetSize(GetNDof());
4850 for (int p = 0; p < GetNP(); p++)
4852 p2g.SetPatchDofMap(p, kv);
4853 NURBSPatch &patch = *patches[p];
4854 MFEM_ASSERT(vdim+1 == patch.GetNC(), "");
4856 for (int k = 0; k < kv[2]->GetNCP(); k++)
4858 for (int j = 0; j < kv[1]->GetNCP(); j++)
4860 for (int i = 0; i < kv[0]->GetNCP(); i++)
4862 const int l = p2g(i,j,k);
4863 for (int d = 0; d < vdim; d++)
4865 coords(l*vdim + d) = patch(i,j,k,d)/patch(i,j,k,vdim);
4867 weights(l) = patch(i,j,k,vdim);
4875void NURBSExtension::GetElementIJK(int elem, Array<int> & ijk)
4877 MFEM_VERIFY(ijk.Size() == el_to_IJK.NumCols(), "");
4878 el_to_IJK.GetRow(elem, ijk);
4881void NURBSExtension::SetPatchToElements()
4883 const int np = GetNP();
4884 patch_to_el.resize(np);
4886 for (int e=0; e<el_to_patch.Size(); ++e)
4888 patch_to_el[el_to_patch[e]].Append(e);
4892void NURBSExtension::SetPatchToBdrElements()
4894 const int nbp = GetNBP();
4895 patch_to_bel.resize(nbp);
4897 for (int e=0; e<bel_to_patch.Size(); ++e)
4899 patch_to_bel[bel_to_patch[e]].Append(e);
4903const Array<int>& NURBSExtension::GetPatchElements(int patch)
4905 MFEM_ASSERT(patch_to_el.size() > 0, "patch_to_el not set
");
4907 return patch_to_el[patch];
4910const Array<int>& NURBSExtension::GetPatchBdrElements(int patch)
4912 MFEM_ASSERT(patch_to_bel.size() > 0, "patch_to_el not set
");
4914 return patch_to_bel[patch];
4917NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim_,
4918 const real_t* control_points)
4921 kv[0] = new KnotVector(*kv0);
4922 kv[1] = new KnotVector(*kv1);
4924 memcpy(data, control_points, sizeof (real_t) * ni * nj * dim_);
4927NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
4928 const KnotVector *kv2, int dim_,
4929 const real_t* control_points)
4932 kv[0] = new KnotVector(*kv0);
4933 kv[1] = new KnotVector(*kv1);
4934 kv[2] = new KnotVector(*kv2);
4936 memcpy(data, control_points, sizeof (real_t) * ni * nj * nk * dim_);
4939NURBSPatch::NURBSPatch(Array<const KnotVector *> &kv_, int dim_,
4940 const real_t* control_points)
4942 kv.SetSize(kv_.Size());
4944 for (int i = 0; i < kv.Size(); i++)
4946 kv[i] = new KnotVector(*kv_[i]);
4947 n *= kv[i]->GetNCP();
4950 memcpy(data, control_points, sizeof(real_t)*n);
4954ParNURBSExtension::ParNURBSExtension(const ParNURBSExtension &orig)
4955 : NURBSExtension(orig),
4956 partitioning(orig.partitioning),
4958 ldof_group(orig.ldof_group)
4962ParNURBSExtension::ParNURBSExtension(MPI_Comm comm, NURBSExtension *parent,
4963 const int *partitioning_,
4964 const Array<bool> &active_bel)
4967 if (parent->NumOfActiveElems < parent->NumOfElements)
4969 // SetActive (BuildGroups?) and the way the weights are copied
4970 // do not support this case
4972 " all elements in the parent must be active!
");
4975 patchTopo = parent->patchTopo;
4976 // steal ownership of patchTopo from the 'parent' NURBS extension
4977 if (!parent->own_topo)
4980 " parent does not own the patch topology!
");
4983 parent->own_topo = false;
4985 parent->edge_to_knot.Copy(edge_to_knot);
4987 parent->GetOrders().Copy(mOrders);
4988 mOrder = parent->GetOrder();
4990 NumOfKnotVectors = parent->GetNKV();
4991 knotVectors.SetSize(NumOfKnotVectors);
4992 for (int i = 0; i < NumOfKnotVectors; i++)
4994 knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
4996 CreateComprehensiveKV();
5002 // copy 'partitioning_' to 'partitioning'
5003 partitioning.SetSize(GetGNE());
5004 for (int i = 0; i < GetGNE(); i++)
5006 partitioning[i] = partitioning_[i];
5008 SetActive(partitioning, active_bel);
5010 GenerateActiveVertices();
5011 GenerateElementDofTable();
5012 // GenerateActiveBdrElems(); // done by SetActive for now
5013 GenerateBdrElementDofTable();
5015 Table *serial_elem_dof = parent->GetElementDofTable();
5016 BuildGroups(partitioning, *serial_elem_dof);
5018 weights.SetSize(GetNDof());
5019 // copy weights from parent
5020 for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
5022 if (activeElem[gel])
5024 int ndofs = el_dof->RowSize(lel);
5025 int *ldofs = el_dof->GetRow(lel);
5026 int *gdofs = serial_elem_dof->GetRow(gel);
5027 for (int i = 0; i < ndofs; i++)
5029 weights(ldofs[i]) = parent->weights(gdofs[i]);
5036ParNURBSExtension::ParNURBSExtension(NURBSExtension *parent,
5037 const ParNURBSExtension *par_parent)
5038 : gtopo(par_parent->gtopo.GetComm())
5040 // steal all data from parent
5041 mOrder = parent->mOrder;
5042 Swap(mOrders, parent->mOrders);
5044 patchTopo = parent->patchTopo;
5045 own_topo = parent->own_topo;
5046 parent->own_topo = false;
5048 Swap(edge_to_knot, parent->edge_to_knot);
5050 NumOfKnotVectors = parent->NumOfKnotVectors;
5051 Swap(knotVectors, parent->knotVectors);
5052 Swap(knotVectorsCompr, parent->knotVectorsCompr);
5054 NumOfVertices = parent->NumOfVertices;
5055 NumOfElements = parent->NumOfElements;
5056 NumOfBdrElements = parent->NumOfBdrElements;
5057 NumOfDofs = parent->NumOfDofs;
5059 Swap(v_meshOffsets, parent->v_meshOffsets);
5060 Swap(e_meshOffsets, parent->e_meshOffsets);
5061 Swap(f_meshOffsets, parent->f_meshOffsets);
5062 Swap(p_meshOffsets, parent->p_meshOffsets);
5064 Swap(v_spaceOffsets, parent->v_spaceOffsets);
5065 Swap(e_spaceOffsets, parent->e_spaceOffsets);
5066 Swap(f_spaceOffsets, parent->f_spaceOffsets);
5067 Swap(p_spaceOffsets, parent->p_spaceOffsets);
5069 Swap(d_to_d, parent->d_to_d);
5070 Swap(master, parent->master);
5071 Swap(slave, parent->slave);
5073 NumOfActiveVertices = parent->NumOfActiveVertices;
5074 NumOfActiveElems = parent->NumOfActiveElems;
5075 NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
5076 NumOfActiveDofs = parent->NumOfActiveDofs;
5078 Swap(activeVert, parent->activeVert);
5079 Swap(activeElem, parent->activeElem);
5080 Swap(activeBdrElem, parent->activeBdrElem);
5081 Swap(activeDof, parent->activeDof);
5083 el_dof = parent->el_dof;
5084 bel_dof = parent->bel_dof;
5085 parent->el_dof = parent->bel_dof = NULL;
5087 Swap(el_to_patch, parent->el_to_patch);
5088 Swap(bel_to_patch, parent->bel_to_patch);
5089 Swap(el_to_IJK, parent->el_to_IJK);
5090 Swap(bel_to_IJK, parent->bel_to_IJK);
5092 Swap(weights, parent->weights);
5093 MFEM_VERIFY(!parent->HavePatches(), "");
5097 MFEM_VERIFY(par_parent->partitioning,
5100 // Support for the case when 'parent' is not a local NURBSExtension, i.e.
5101 // NumOfActiveElems is not the same as in 'par_parent'. In that case, we
5102 // assume 'parent' is a global NURBSExtension, i.e. all elements are active.
5103 bool extract_weights = false;
5104 if (NumOfActiveElems != par_parent->NumOfActiveElems)
5106 MFEM_ASSERT(NumOfActiveElems == NumOfElements, "internal error
");
5108 SetActive(par_parent->partitioning, par_parent->activeBdrElem);
5109 GenerateActiveVertices();
5111 el_to_patch.DeleteAll();
5112 el_to_IJK.DeleteAll();
5113 GenerateElementDofTable();
5114 // GenerateActiveBdrElems(); // done by SetActive for now
5116 bel_to_patch.DeleteAll();
5117 bel_to_IJK.DeleteAll();
5118 GenerateBdrElementDofTable();
5119 extract_weights = true;
5122 Table *glob_elem_dof = GetGlobalElementDofTable();
5123 BuildGroups(par_parent->partitioning, *glob_elem_dof);
5124 if (extract_weights)
5126 Vector glob_weights;
5127 Swap(weights, glob_weights);
5128 weights.SetSize(GetNDof());
5129 // Copy the local 'weights' from the 'glob_weights'.
5130 // Assumption: the local element ids follow the global ordering.
5131 for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
5133 if (activeElem[gel])
5135 int ndofs = el_dof->RowSize(lel);
5136 int *ldofs = el_dof->GetRow(lel);
5137 int *gdofs = glob_elem_dof->GetRow(gel);
5138 for (int i = 0; i < ndofs; i++)
5140 weights(ldofs[i]) = glob_weights(gdofs[i]);
5146 delete glob_elem_dof;
5149Table *ParNURBSExtension::GetGlobalElementDofTable()
5151 if (Dimension() == 1)
5153 return Get1DGlobalElementDofTable();
5155 else if (Dimension() == 2)
5157 return Get2DGlobalElementDofTable();
5161 return Get3DGlobalElementDofTable();
5165Table *ParNURBSExtension::Get1DGlobalElementDofTable()
5168 const KnotVector *kv[1];
5169 NURBSPatchMap p2g(this);
5170 Array<Connection> gel_dof_list;
5172 for (int p = 0; p < GetNP(); p++)
5174 p2g.SetPatchDofMap(p, kv);
5177 const int ord0 = kv[0]->GetOrder();
5179 for (int i = 0; i < kv[0]->GetNKS(); i++)
5181 if (kv[0]->isElement(i))
5183 Connection conn(el,0);
5184 for (int ii = 0; ii <= ord0; ii++)
5186 conn.to = DofMap(p2g(i+ii));
5187 gel_dof_list.Append(conn);
5193 // We must NOT sort gel_dof_list in this case.
5194 return (new Table(GetGNE(), gel_dof_list));
5197Table *ParNURBSExtension::Get2DGlobalElementDofTable()
5200 const KnotVector *kv[2];
5201 NURBSPatchMap p2g(this);
5202 Array<Connection> gel_dof_list;
5204 for (int p = 0; p < GetNP(); p++)
5206 p2g.SetPatchDofMap(p, kv);
5209 const int ord0 = kv[0]->GetOrder();
5210 const int ord1 = kv[1]->GetOrder();
5211 for (int j = 0; j < kv[1]->GetNKS(); j++)
5213 if (kv[1]->isElement(j))
5215 for (int i = 0; i < kv[0]->GetNKS(); i++)
5217 if (kv[0]->isElement(i))
5219 Connection conn(el,0);
5220 for (int jj = 0; jj <= ord1; jj++)
5222 for (int ii = 0; ii <= ord0; ii++)
5224 conn.to = DofMap(p2g(i+ii,j+jj));
5225 gel_dof_list.Append(conn);
5234 // We must NOT sort gel_dof_list in this case.
5235 return (new Table(GetGNE(), gel_dof_list));
5238Table *ParNURBSExtension::Get3DGlobalElementDofTable()
5241 const KnotVector *kv[3];
5242 NURBSPatchMap p2g(this);
5243 Array<Connection> gel_dof_list;
5245 for (int p = 0; p < GetNP(); p++)
5247 p2g.SetPatchDofMap(p, kv);
5250 const int ord0 = kv[0]->GetOrder();
5251 const int ord1 = kv[1]->GetOrder();
5252 const int ord2 = kv[2]->GetOrder();
5253 for (int k = 0; k < kv[2]->GetNKS(); k++)
5255 if (kv[2]->isElement(k))
5257 for (int j = 0; j < kv[1]->GetNKS(); j++)
5259 if (kv[1]->isElement(j))
5261 for (int i = 0; i < kv[0]->GetNKS(); i++)
5263 if (kv[0]->isElement(i))
5265 Connection conn(el,0);
5266 for (int kk = 0; kk <= ord2; kk++)
5268 for (int jj = 0; jj <= ord1; jj++)
5270 for (int ii = 0; ii <= ord0; ii++)
5272 conn.to = DofMap(p2g(i+ii,j+jj,k+kk));
5273 gel_dof_list.Append(conn);
5285 // We must NOT sort gel_dof_list in this case.
5286 return (new Table(GetGNE(), gel_dof_list));
5289void ParNURBSExtension::SetActive(const int *partition,
5290 const Array<bool> &active_bel)
5292 activeElem.SetSize(GetGNE());
5294 NumOfActiveElems = 0;
5295 const int MyRank = gtopo.MyRank();
5296 for (int i = 0; i < GetGNE(); i++)
5297 if (partition[i] == MyRank)
5299 activeElem[i] = true;
5303 active_bel.Copy(activeBdrElem);
5304 NumOfActiveBdrElems = 0;
5305 for (int i = 0; i < GetGNBE(); i++)
5306 if (activeBdrElem[i])
5308 NumOfActiveBdrElems++;
5312void ParNURBSExtension::BuildGroups(const int *partition,
5313 const Table &elem_dof)
5317 ListOfIntegerSets groups;
5320 Transpose(elem_dof, dof_proc); // dof_proc is dof_elem
5321 // convert elements to processors
5322 for (int i = 0; i < dof_proc.Size_of_connections(); i++)
5324 dof_proc.GetJ()[i] = partition[dof_proc.GetJ()[i]];
5327 // the first group is the local one
5328 int MyRank = gtopo.MyRank();
5329 group.Recreate(1, &MyRank);
5330 groups.Insert(group);
5333 ldof_group.SetSize(GetNDof());
5334 for (int d = 0; d < GetNTotalDof(); d++)
5337 group.Recreate(dof_proc.RowSize(d), dof_proc.GetRow(d));
5338 ldof_group[dof] = groups.Insert(group);
5343 gtopo.Create(groups, 1822);
5345#endif // MFEM_USE_MPI
5348void NURBSPatchMap::GetPatchKnotVectors(int p, const KnotVector *kv[])
5350 Ext->patchTopo->GetElementVertices(p, verts);
5352 if (Ext->Dimension() == 1)
5354 kv[0] = Ext->knotVectorsCompr[Ext->Dimension()*p];
5356 else if (Ext->Dimension() == 2)
5358 Ext->patchTopo->GetElementEdges(p, edges, oedge);
5360 kv[0] = Ext->knotVectorsCompr[Ext->Dimension()*p];
5361 kv[1] = Ext->knotVectorsCompr[Ext->Dimension()*p + 1];
5363 else if (Ext->Dimension() == 3)
5365 Ext->patchTopo->GetElementEdges(p, edges, oedge);
5366 Ext->patchTopo->GetElementFaces(p, faces, oface);
5368 kv[0] = Ext->knotVectorsCompr[Ext->Dimension()*p];
5369 kv[1] = Ext->knotVectorsCompr[Ext->Dimension()*p + 1];
5370 kv[2] = Ext->knotVectorsCompr[Ext->Dimension()*p + 2];
5375void NURBSPatchMap::GetBdrPatchKnotVectors(int p, const KnotVector *kv[],
5378 Ext->patchTopo->GetBdrElementVertices(p, verts);
5380 if (Ext->Dimension() == 2)
5382 Ext->patchTopo->GetBdrElementEdges(p, edges, oedge);
5383 kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]);
5386 else if (Ext->Dimension() == 3)
5389 Ext->patchTopo->GetBdrElementEdges(p, edges, oedge);
5390 Ext->patchTopo->GetBdrElementFace(p, &faces[0], &opatch);
5392 kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]);
5393 kv[1] = Ext->KnotVec(edges[1], oedge[1], &okv[1]);
5398void NURBSPatchMap::SetPatchVertexMap(int p, const KnotVector *kv[])
5400 GetPatchKnotVectors(p, kv);
5402 I = kv[0]->GetNE() - 1;
5404 for (int i = 0; i < verts.Size(); i++)
5406 verts[i] = Ext->v_meshOffsets[verts[i]];
5409 if (Ext->Dimension() >= 2)
5411 J = kv[1]->GetNE() - 1;
5412 for (int i = 0; i < edges.Size(); i++)
5414 edges[i] = Ext->e_meshOffsets[edges[i]];
5417 if (Ext->Dimension() == 3)
5419 K = kv[2]->GetNE() - 1;
5421 for (int i = 0; i < faces.Size(); i++)
5423 faces[i] = Ext->f_meshOffsets[faces[i]];
5427 pOffset = Ext->p_meshOffsets[p];
5430void NURBSPatchMap::SetPatchDofMap(int p, const KnotVector *kv[])
5432 GetPatchKnotVectors(p, kv);
5434 I = kv[0]->GetNCP() - 2;
5436 for (int i = 0; i < verts.Size(); i++)
5438 verts[i] = Ext->v_spaceOffsets[verts[i]];
5440 if (Ext->Dimension() >= 2)
5442 J = kv[1]->GetNCP() - 2;
5443 for (int i = 0; i < edges.Size(); i++)
5445 edges[i] = Ext->e_spaceOffsets[edges[i]];
5448 if (Ext->Dimension() == 3)
5450 K = kv[2]->GetNCP() - 2;
5452 for (int i = 0; i < faces.Size(); i++)
5454 faces[i] = Ext->f_spaceOffsets[faces[i]];
5458 pOffset = Ext->p_spaceOffsets[p];
5461void NURBSPatchMap::SetBdrPatchVertexMap(int p, const KnotVector *kv[],
5464 GetBdrPatchKnotVectors(p, kv, okv);
5466 for (int i = 0; i < verts.Size(); i++)
5468 verts[i] = Ext->v_meshOffsets[verts[i]];
5471 if (Ext->Dimension() == 1)
5475 else if (Ext->Dimension() == 2)
5477 I = kv[0]->GetNE() - 1;
5478 pOffset = Ext->e_meshOffsets[edges[0]];
5480 else if (Ext->Dimension() == 3)
5482 I = kv[0]->GetNE() - 1;
5483 J = kv[1]->GetNE() - 1;
5485 for (int i = 0; i < edges.Size(); i++)
5487 edges[i] = Ext->e_meshOffsets[edges[i]];
5490 pOffset = Ext->f_meshOffsets[faces[0]];
5494void NURBSPatchMap::SetBdrPatchDofMap(int p, const KnotVector *kv[], int *okv)
5496 GetBdrPatchKnotVectors(p, kv, okv);
5498 for (int i = 0; i < verts.Size(); i++)
5500 verts[i] = Ext->v_spaceOffsets[verts[i]];
5503 if (Ext->Dimension() == 1)
5507 else if (Ext->Dimension() == 2)
5509 I = kv[0]->GetNCP() - 2;
5510 pOffset = Ext->e_spaceOffsets[edges[0]];
5512 else if (Ext->Dimension() == 3)
5514 I = kv[0]->GetNCP() - 2;
5515 J = kv[1]->GetNCP() - 2;
5517 for (int i = 0; i < edges.Size(); i++)
5519 edges[i] = Ext->e_spaceOffsets[edges[i]];
5522 pOffset = Ext->f_spaceOffsets[faces[0]];
int Size() const
Return the logical size of the array.
T Sum() const
Return the sum of all the array entries using the '+'' operator for class 'T'.
A vector of knots in one dimension, with B-spline basis functions of a prescribed order.
std::shared_ptr< SpacingFunction > spacing
Function to define the distribution of knots for any number of knot spans.
void PrintFunctions(std::ostream &os, int samples=11) const
Prints the non-zero shape functions and their first and second derivatives associated with the KnotVe...
int NumOfElements
Number of elements, defined by distinct knots.
void CalcDnShape(Vector &gradn, int n, int i, real_t xi) const
Calculate n-th derivatives (order n) of the nonvanishing shape function values in grad for the elemen...
int Order
Order of the B-spline basis functions.
KnotVector & operator=(const KnotVector &kv)
void UniformRefinement(Vector &newknots, int rf) const
Uniformly refine by factor rf, by inserting knots in each span.
bool isElement(int i) const
Return whether the knot index Order plus i is the beginning of an element.
void CalcShape(Vector &shape, int i, real_t xi) const
Calculate the nonvanishing shape function values in shape for the element corresponding to knot index...
Vector knot
Stores the values of all knots.
int GetNKS() const
Return the number of control points minus the order. This is not the number of knot spans,...
void GetElements()
Count the number of elements.
KnotVector * DegreeElevate(int t) const
Return a new KnotVector with elevated degree by repeating the endpoints of the KnotVector.
void Refinement(Vector &newknots, int rf) const
Refine with refinement factor rf.
KnotVector()
Create an empty KnotVector.
static const int MaxOrder
void CalcD2Shape(Vector &grad2, int i, real_t xi) const
Calculate second-order shape function derivatives, using CalcDnShape.
int GetNCP() const
Return the number of control points.
int NumOfControlPoints
Number of control points.
void CalcDShape(Vector &grad, int i, real_t xi) const
Calculate derivatives of the nonvanishing shape function values in grad for the element corresponding...
int Size() const
Return the number of knots, including multiplicities.
void Flip()
Reverse the knots.
void Difference(const KnotVector &kv, Vector &diff) const
int GetCoarseningFactor() const
Vector GetFineKnots(const int cf) const
int GetNE() const
Return the number of elements, defined by distinct knots.
void Print(std::ostream &os) const
Print the order, number of control points, and knots.
NURBSExtension generally contains multiple NURBSPatch objects spanning an entire Mesh....
A NURBS patch can be 1D, 2D, or 3D, and is defined as a tensor product of KnotVectors.
Parallel version of NURBSExtension.
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
int index(int i, int j, int nx, int ny)
void mfem_error(const char *msg)
NURBSPatch * Interpolate(NURBSPatch &p1, NURBSPatch &p2)
NURBSPatch * Revolve3D(NURBSPatch &patch, real_t n[], real_t ang, int times)