12 #include "../fem/fem.hpp"
14 #if defined(_MSC_VER) && (_MSC_VER < 1800)
16 #define copysign _copysign
28 input >> Order >> NumOfControlPoints;
29 knot.Load(input, NumOfControlPoints + Order + 1);
36 NumOfControlPoints = NCP;
37 knot.SetSize(NumOfControlPoints + Order + 1);
58 " Parent KnotVector order higher than child");
60 int nOrder = Order + t;
63 for (
int i = 0; i <= nOrder; i++)
65 (*newkv)[i] = knot(0);
67 for (
int i = nOrder + 1; i < newkv->
GetNCP(); i++)
69 (*newkv)[i] = knot(i - t);
71 for (
int i = 0; i <= nOrder; i++)
73 (*newkv)[newkv->
GetNCP() + i] = knot(knot.Size()-1);
83 newknots.
SetSize(NumOfElements);
85 for (
int i = 0; i < knot.Size()-1; i++)
87 if (knot(i) != knot(i+1))
89 newknots(j) = 0.5*(knot(i) + knot(i+1));
98 for (
int i = Order; i < NumOfControlPoints; i++)
100 if (knot(i) != knot(i+1))
109 double apb = knot(0) + knot(knot.Size()-1);
111 int ns = (NumOfControlPoints - Order)/2;
112 for (
int i = 1; i <= ns; i++)
114 double tmp = apb - knot(Order + i);
115 knot(Order + i) = apb - knot(NumOfControlPoints - i);
116 knot(NumOfControlPoints - i) = tmp;
122 out << Order <<
' ' << NumOfControlPoints <<
' ';
123 knot.Print(out, knot.Size());
130 int ip = (i >= 0) ? (i + p) : (-1 - i + p);
131 double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), saved, tmp;
132 double left[MaxOrder+1], right[MaxOrder+1];
137 mfem_error(
"KnotVector::CalcShape : Order > MaxOrder!");
142 for (
int j = 1; j <= p; ++j)
144 left[j] = u - knot(ip+1-j);
145 right[j] = knot(ip+j) - u;
147 for (
int r = 0; r < j; ++r)
149 tmp = shape(r)/(right[r+1] + left[j-r]);
150 shape(r) = saved + right[r+1]*tmp;
151 saved = left[j-r]*tmp;
160 int p = Order, rk, pk;
161 int ip = (i >= 0) ? (i + p) : (-1 - i + p);
162 double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), temp, saved, d;
163 double ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1], right[MaxOrder+1];
168 mfem_error(
"KnotVector::CalcDShape : Order > MaxOrder!");
173 for (
int j = 1; j <= p; j++)
175 left[j] = u - knot(ip-j+1);
176 right[j] = knot(ip+j) - u;
178 for (
int r = 0; r < j; r++)
180 ndu[j][r] = right[r+1] + left[j-r];
181 temp = ndu[r][j-1]/ndu[j][r];
182 ndu[r][j] = saved + right[r+1]*temp;
183 saved = left[j-r]*temp;
188 for (
int r = 0; r <= p; ++r)
195 d = ndu[rk][pk]/ndu[p][rk];
199 d -= ndu[r][pk]/ndu[p][r];
206 grad *= p*(knot(ip+1) - knot(ip));
210 grad *= p*(knot(ip) - knot(ip+1));
218 if (u == knot(NumOfControlPoints+Order))
220 mid = NumOfControlPoints;
225 high = NumOfControlPoints + 1;
226 mid = (low + high)/2;
227 while ( (u < knot(mid-1)) || (u > knot(mid)) )
237 mid = (low + high)/2;
248 " Can not compare knot vectors with different orders!");
251 int s = kv.
Size() - Size();
262 for (
int j = 0; j < kv.
Size(); j++)
264 if (knot(i) == kv[j])
284 ni = kv[0]->GetNCP();
285 nj = kv[1]->GetNCP();
288 data =
new double[ni*nj*Dim];
291 for (
int i = 0; i < ni*nj*Dim; i++)
297 else if (kv.Size() == 3)
299 ni = kv[0]->GetNCP();
300 nj = kv[1]->GetNCP();
301 nk = kv[2]->GetNCP();
303 data =
new double[ni*nj*nk*Dim];
306 for (
int i = 0; i < ni*nj*nk*Dim; i++)
314 mfem_error(
"NURBSPatch::init : Wrond dimension of knotvectors!");
320 int pdim,
dim, size = 1;
323 input >> ws >> ident >> pdim;
325 for (
int i = 0; i < pdim; i++)
328 size *= kv[i]->GetNCP();
331 input >> ws >> ident >>
dim;
334 input >> ws >> ident;
335 if (ident ==
"controlpoints" || ident ==
"controlpoints_homogeneous")
337 for (
int j = 0, i = 0; i < size; i++)
338 for (
int d = 0; d <=
dim; d++, j++)
345 for (
int j = 0, i = 0; i < size; i++)
347 for (
int d = 0; d <=
dim; d++)
351 for (
int d = 0; d <
dim; d++)
353 data[j+d] *= data[j+
dim];
380 kv.SetSize(kv_.
Size());
381 for (
int i = 0; i < kv.Size(); i++)
390 kv.SetSize(parent->
kv.Size());
391 for (
int i = 0; i < kv.Size(); i++)
410 for (
int i = 0; i < kv.Size(); i++)
412 if (kv[i]) {
delete kv[i]; }
436 for (
int i = 0; i < kv.Size(); i++)
438 if (kv[i]) {
delete kv[i]; }
446 out <<
"knotvectors\n" << kv.Size() <<
'\n';
447 for (
int i = 0; i < kv.Size(); i++)
450 size *= kv[i]->GetNCP();
453 out <<
"\ndimension\n" << Dim - 1
454 <<
"\n\ncontrolpoints\n";
455 for (
int j = 0, i = 0; i < size; i++)
458 for (
int d = 1; d < Dim; d++)
460 out <<
' ' << data[j++];
486 cerr <<
"NURBSPatch::SetLoopDirection :\n"
487 " Direction error in 2D patch, dir = " << dir <<
'\n';
516 cerr <<
"NURBSPatch::SetLoopDirection :\n"
517 " Direction error in 3D patch, dir = " << dir <<
'\n';
528 for (
int dir = 0; dir < kv.Size(); dir++)
530 kv[dir]->UniformRefinement(newknots);
531 KnotInsert(dir, newknots);
537 for (
int dir = 0; dir < kv.Size(); dir++)
539 KnotInsert(dir, *newkv[dir]);
545 if (dir >= kv.Size() || dir < 0)
547 mfem_error(
"NURBSPatch::KnotInsert : Incorrect direction!");
550 int t = newkv.
GetOrder() - kv[dir]->GetOrder();
554 DegreeElevate(dir, t);
558 mfem_error(
"NURBSPatch::KnotInsert : Incorrect order!");
562 GetKV(dir)->Difference(newkv, diff);
565 KnotInsert(dir, diff);
572 if (dir >= kv.Size() || dir < 0)
574 mfem_error(
"NURBSPatch::KnotInsert : Incorrect direction!");
586 if (size != newp.SetLoopDirection(dir))
588 mfem_error(
"NURBSPatch::KnotInsert : Size mismatch!");
591 int rr = knot.
Size() - 1;
597 for (
int j = 0; j <= a; j++)
601 for (
int j = b+pl; j <= ml+pl; j++)
603 newkv[j+rr+1] = oldkv[j];
605 for (
int k = 0; k <= (a-pl); k++)
607 for (
int ll = 0; ll < size; ll++)
609 newp(k,ll) = oldp(k,ll);
612 for (
int k = (b-1); k < ml; k++)
614 for (
int ll = 0; ll < size; ll++)
616 newp(k+rr+1,ll) = oldp(k,ll);
623 for (
int j = rr; j >= 0; j--)
625 while ( (knot(j) <= oldkv[i]) && (i > a) )
628 for (
int ll = 0; ll < size; ll++)
630 newp(k-pl-1,ll) = oldp(i-pl-1,ll);
637 for (
int ll = 0; ll < size; ll++)
639 newp(k-pl-1,ll) = newp(k-pl,ll);
642 for (
int l = 1; l <= pl; l++)
645 double alfa = newkv[k+l] - knot(j);
646 if (fabs(alfa) == 0.0)
648 for (
int ll = 0; ll < size; ll++)
650 newp(ind-1,ll) = newp(ind,ll);
655 alfa = alfa/(newkv[k+l] - oldkv[i-pl+l]);
656 for (
int ll = 0; ll < size; ll++)
658 newp(ind-1,ll) = alfa*newp(ind-1,ll) + (1.0-alfa)*newp(ind,ll);
674 for (
int dir = 0; dir < kv.Size(); dir++)
676 DegreeElevate(dir, t);
683 if (dir >= kv.Size() || dir < 0)
685 mfem_error(
"NURBSPatch::DegreeElevate : Incorrect direction!");
688 int i, j, k, kj, mpi, mul, mh, kind, cind, first, last;
689 int r, a, b, oldr, save, s, tr, lbz, rbz, l;
690 double inv, ua, ub, numer, alf, den, bet, gam;
701 if (size != newp.SetLoopDirection(dir))
703 mfem_error(
"NURBSPatch::DegreeElevate : Size mismatch!");
721 for (i = 0; i <= ph; i++)
723 binom(i,0) = binom(i,i) = 1;
724 for (j = 1; j < i; j++)
726 binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
733 for (i = 1; i <= ph2; i++)
735 inv = 1.0/binom(ph,i);
737 for (j = max(0,i-t); j <= mpi; j++)
739 bezalfs(i,j) = inv*binom(p,j)*binom(t,i-j);
744 for (i = ph2+1; i < ph; i++)
747 for (j = max(0,i-t); j <= mpi; j++)
749 bezalfs(i,j) = bezalfs(ph-i,p-j);
760 for (l = 0; l < size; l++)
762 newp(0,l) = oldp(0,l);
764 for (i = 0; i <= ph; i++)
769 for (i = 0; i <= p; i++)
771 for (l = 0; l < size; l++)
773 bpts(i,l) = oldp(i,l);
780 while (b < m && oldkv[b] == oldkv[b+1]) { b++; }
788 if (oldr > 0) { lbz = (oldr+2)/2; }
791 if (r > 0) { rbz = ph-(r+1)/2; }
797 for (k = p ; k > mul; k--)
799 alphas[k-mul-1] = numer/(oldkv[a+k]-ua);
802 for (j = 1; j <= r; j++)
806 for (k = p; k >= s; k--)
808 for (l = 0; l < size; l++)
809 bpts(k,l) = (alphas[k-s]*bpts(k,l) +
810 (1.0-alphas[k-s])*bpts(k-1,l));
812 for (l = 0; l < size; l++)
814 nextbpts(save,l) = bpts(p,l);
819 for (i = lbz; i <= ph; i++)
821 for (l = 0; l < size; l++)
826 for (j = max(0,i-t); j <= mpi; j++)
828 for (l = 0; l < size; l++)
830 ebpts(i,l) += bezalfs(i,j)*bpts(j,l);
840 bet = (ub-newkv[kind-1])/den;
842 for (tr = 1; tr < oldr; tr++)
851 alf = (ub-newkv[i])/(ua-newkv[i]);
852 for (l = 0; l < size; l++)
854 newp(i,l) = alf*newp(i,l)-(1.0-alf)*newp(i-1,l);
859 if ((j-tr) <= (kind-ph+oldr))
861 gam = (ub-newkv[j-tr])/den;
862 for (l = 0; l < size; l++)
864 ebpts(kj,l) = gam*ebpts(kj,l) + (1.0-gam)*ebpts(kj+1,l);
869 for (l = 0; l < size; l++)
871 ebpts(kj,l) = bet*ebpts(kj,l) + (1.0-bet)*ebpts(kj+1,l);
886 for (i = 0; i < (ph-oldr); i++)
892 for (j = lbz; j <= rbz; j++)
894 for (l = 0; l < size; l++)
896 newp(cind,l) = ebpts(j,l);
903 for (j = 0; j <r; j++)
904 for (l = 0; l < size; l++)
906 bpts(j,l) = nextbpts(j,l);
909 for (j = r; j <= p; j++)
910 for (l = 0; l < size; l++)
912 bpts(j,l) = oldp(b-p+j,l);
921 for (i = 0; i <= ph; i++)
934 int size = SetLoopDirection(dir);
936 for (
int id = 0;
id < nd/2;
id++)
937 for (
int i = 0; i < size; i++)
939 Swap<double>((*this)(id,i), (*
this)(nd-1-id,i));
946 if (abs(dir1-dir2) == 2)
948 " directions 0 and 2 are not supported!");
953 Swap<KnotVector *>(nkv[dir1], nkv[dir2]);
956 int size = SetLoopDirection(dir1);
957 newpatch->SetLoopDirection(dir2);
959 for (
int id = 0;
id < nd;
id++)
960 for (
int i = 0; i < size; i++)
962 (*newpatch)(id,i) = (*
this)(id,i);
972 double l2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
975 if (fabs(angle) == M_PI_2)
977 s = r*copysign(1., angle);
981 else if (fabs(angle) == M_PI)
996 T(0,0) = (n[0]*n[0] + (n[1]*n[1] + n[2]*n[2])*c)/l2;
997 T(0,1) = -(n[0]*n[1]*c1)/l2 - (n[2]*s)/l;
998 T(0,2) = -(n[0]*n[2]*c1)/l2 + (n[1]*s)/l;
999 T(1,0) = -(n[0]*n[1]*c1)/l2 + (n[2]*s)/l;
1000 T(1,1) = (n[1]*n[1] + (n[0]*n[0] + n[2]*n[2])*c)/l2;
1001 T(1,2) = -(n[1]*n[2]*c1)/l2 - (n[0]*s)/l;
1002 T(2,0) = -(n[0]*n[2]*c1)/l2 - (n[1]*s)/l;
1003 T(2,1) = -(n[1]*n[2]*c1)/l2 + (n[0]*s)/l;
1004 T(2,2) = (n[2]*n[2] + (n[0]*n[0] + n[1]*n[1])*c)/l2;
1011 mfem_error(
"NURBSPatch::Rotate3D : not a NURBSPatch in 3D!");
1017 Get3DRotationMatrix(n, angle, 1., T);
1020 for (
int i = 0; i < kv.Size(); i++)
1022 size *= kv[i]->GetNCP();
1025 for (
int i = 0; i < size; i++)
1037 for (
int dir = 0; dir < kv.Size(); dir++)
1038 if (maxd < kv[dir]->GetOrder())
1040 maxd = kv[dir]->GetOrder();
1043 for (
int dir = 0; dir < kv.Size(); dir++)
1044 if (maxd > kv[dir]->GetOrder())
1046 DegreeElevate(dir, maxd - kv[dir]->GetOrder());
1054 if (p1.
kv.Size() != p2.
kv.Size() || p1.
Dim != p2.
Dim)
1056 mfem_error(
"Interpolate(NURBSPatch &, NURBSPatch &)");
1059 int size = 1,
dim = p1.
Dim;
1062 for (
int i = 0; i < p1.
kv.Size(); i++)
1064 if (p1.
kv[i]->GetOrder() < p2.
kv[i]->GetOrder())
1075 size *= kv[i]->GetNCP();
1080 nkv[0] = nkv[1] = 0.0;
1081 nkv[2] = nkv[3] = 1.0;
1087 for (
int i = 0; i < size; i++)
1089 for (
int d = 0; d <
dim; d++)
1091 patch->
data[i*dim+d] = p1.
data[i*dim+d];
1092 patch->
data[(i+size)*dim+d] = p2.
data[i*dim+d];
1103 mfem_error(
"Revolve3D(NURBSPatch &, double [], double)");
1109 for (
int i = 0; i < patch.
kv.Size(); i++)
1111 nkv[i] = patch.
kv[i];
1112 size *= nkv[i]->GetNCP();
1117 lkv[0] = lkv[1] = lkv[2] = 0.0;
1118 for (
int i = 1; i < times; i++)
1120 lkv[2*i+1] = lkv[2*i+2] = i;
1122 lkv[ns] = lkv[ns+1] = lkv[ns+2] = times;
1128 Vector u(NULL, 3), v(NULL, 3);
1131 double c = cos(ang/2);
1135 double *op = patch.
data, *np;
1136 for (
int i = 0; i < size; i++)
1138 np = newpatch->
data + 4*i;
1139 for (
int j = 0; j < 4; j++)
1143 for (
int j = 0; j < times; j++)
1163 patchTopo =
new Mesh;
1174 input >> ws >> ident;
1175 if (ident ==
"knotvectors")
1177 input >> NumOfKnotVectors;
1178 knotVectors.SetSize(NumOfKnotVectors);
1179 for (
int i = 0; i < NumOfKnotVectors; i++)
1182 if (knotVectors[i]->GetOrder() != knotVectors[0]->GetOrder())
1183 mfem_error(
"NURBSExtension::NURBSExtension :\n"
1184 " Variable orders are not supported!");
1186 Order = knotVectors[0]->GetOrder();
1188 else if (ident ==
"patches")
1190 patches.SetSize(GetNP());
1191 for (
int p = 0; p < patches.Size(); p++)
1197 NumOfKnotVectors = 0;
1198 for (
int i = 0; i < patchTopo->GetNEdges(); i++)
1199 if (NumOfKnotVectors < KnotInd(i))
1201 NumOfKnotVectors = KnotInd(i);
1204 knotVectors.SetSize(NumOfKnotVectors);
1208 for (
int p = 0; p < patches.Size(); p++)
1210 patchTopo->GetElementEdges(p, edges, oedge);
1211 if (Dimension() == 2)
1213 if (knotVectors[KnotInd(edges[0])] == NULL)
1214 knotVectors[KnotInd(edges[0])] =
1216 if (knotVectors[KnotInd(edges[1])] == NULL)
1217 knotVectors[KnotInd(edges[1])] =
1222 if (knotVectors[KnotInd(edges[0])] == NULL)
1223 knotVectors[KnotInd(edges[0])] =
1225 if (knotVectors[KnotInd(edges[3])] == NULL)
1226 knotVectors[KnotInd(edges[3])] =
1228 if (knotVectors[KnotInd(edges[8])] == NULL)
1229 knotVectors[KnotInd(edges[8])] =
1233 Order = knotVectors[0]->GetOrder();
1244 if (patches.Size() == 0)
1246 input >> ws >> ident;
1248 if (patches.Size() == 0 && ident ==
"mesh_elements")
1250 input >> NumOfActiveElems;
1251 activeElem.SetSize(GetGNE());
1254 for (
int i = 0; i < NumOfActiveElems; i++)
1257 activeElem[glob_elem] =
true;
1261 input >> ws >> ident;
1265 NumOfActiveElems = NumOfElements;
1266 activeElem.SetSize(NumOfElements);
1270 GenerateActiveVertices();
1271 GenerateElementDofTable();
1272 GenerateActiveBdrElems();
1273 GenerateBdrElementDofTable();
1275 if (patches.Size() == 0)
1278 if (ident ==
"weights")
1280 weights.Load(input, GetNDof());
1284 weights.SetSize(GetNDof());
1299 NumOfKnotVectors = parent->
GetNKV();
1300 knotVectors.SetSize(NumOfKnotVectors);
1301 for (
int i = 0; i < NumOfKnotVectors; i++)
1320 GenerateElementDofTable();
1321 GenerateBdrElementDofTable();
1323 weights.SetSize(GetNDof());
1332 mfem_error(
"NURBSExtension::NURBSExtension :\n"
1333 " parent does not own the patch topology!");
1342 NumOfKnotVectors = parent->
GetNKV();
1343 knotVectors.SetSize(NumOfKnotVectors);
1344 for (
int i = 0; i < NumOfKnotVectors; i++)
1354 NumOfActiveElems = NumOfElements;
1355 activeElem.SetSize(NumOfElements);
1358 GenerateActiveVertices();
1359 GenerateElementDofTable();
1360 GenerateActiveBdrElems();
1361 GenerateBdrElementDofTable();
1363 weights.SetSize(GetNDof());
1364 MergeWeights(mesh_array, num_pieces);
1369 if (patches.Size() == 0)
1375 for (
int i = 0; i < knotVectors.Size(); i++)
1377 delete knotVectors[i];
1380 for (
int i = 0; i < patches.Size(); i++)
1393 patchTopo->PrintTopo(out, edge_to_knot);
1394 if (patches.Size() == 0)
1396 out <<
"\nknotvectors\n" << NumOfKnotVectors <<
'\n';
1397 for (
int i = 0; i < NumOfKnotVectors; i++)
1399 knotVectors[i]->Print(out);
1402 if (NumOfActiveElems < NumOfElements)
1404 out <<
"\nmesh_elements\n" << NumOfActiveElems <<
'\n';
1405 for (
int i = 0; i < NumOfElements; i++)
1412 out <<
"\nweights\n";
1413 weights.Print(out, 1);
1417 out <<
"\npatches\n";
1418 for (
int p = 0; p < patches.Size(); p++)
1420 out <<
"\n# patch " << p <<
"\n\n";
1421 patches[p]->Print(out);
1429 "NURBS Mesh entity sizes:\n"
1430 "Dimension = " << Dimension() <<
"\n"
1431 "Order = " << GetOrder() <<
"\n"
1432 "NumOfKnotVectors = " << GetNKV() <<
"\n"
1433 "NumOfPatches = " << GetNP() <<
"\n"
1434 "NumOfBdrPatches = " << GetNBP() <<
"\n"
1435 "NumOfVertices = " << GetGNV() <<
"\n"
1436 "NumOfElements = " << GetGNE() <<
"\n"
1437 "NumOfBdrElements = " << GetGNBE() <<
"\n"
1438 "NumOfDofs = " << GetNTotalDof() <<
"\n"
1439 "NumOfActiveVertices = " << GetNV() <<
"\n"
1440 "NumOfActiveElems = " << GetNE() <<
"\n"
1441 "NumOfActiveBdrElems = " << GetNBE() <<
"\n"
1442 "NumOfActiveDofs = " << GetNDof() <<
'\n';
1443 for (
int i = 0; i < NumOfKnotVectors; i++)
1445 out <<
' ' << i + 1 <<
") ";
1446 knotVectors[i]->Print(out);
1453 int vert[8], nv, g_el, nx, ny, nz,
dim = Dimension();
1459 activeVert.SetSize(GetGNV());
1461 for (
int p = 0; p < GetNP(); p++)
1467 nz = (dim == 3) ? p2g.
nz() : 1;
1469 for (
int k = 0; k < nz; k++)
1471 for (
int j = 0; j < ny; j++)
1473 for (
int i = 0; i < nx; i++)
1475 if (activeElem[g_el])
1479 vert[0] = p2g(i, j );
1480 vert[1] = p2g(i+1,j );
1481 vert[2] = p2g(i+1,j+1);
1482 vert[3] = p2g(i, j+1);
1487 vert[0] = p2g(i, j, k);
1488 vert[1] = p2g(i+1,j, k);
1489 vert[2] = p2g(i+1,j+1,k);
1490 vert[3] = p2g(i, j+1,k);
1492 vert[4] = p2g(i, j, k+1);
1493 vert[5] = p2g(i+1,j, k+1);
1494 vert[6] = p2g(i+1,j+1,k+1);
1495 vert[7] = p2g(i, j+1,k+1);
1499 for (
int v = 0; v < nv; v++)
1501 activeVert[vert[v]] = 1;
1510 NumOfActiveVertices = 0;
1511 for (
int i = 0; i < GetGNV(); i++)
1512 if (activeVert[i] == 1)
1514 activeVert[i] = NumOfActiveVertices++;
1520 int dim = Dimension();
1523 activeBdrElem.SetSize(GetGNBE());
1524 if (GetGNE() == GetNE())
1526 activeBdrElem =
true;
1527 NumOfActiveBdrElems = GetGNBE();
1530 activeBdrElem =
false;
1531 NumOfActiveBdrElems = 0;
1544 for (
int i = 0; i < num_pieces; i++)
1550 for (
int lel = 0; lel < lext->
GetNE(); lel++)
1552 int gel = lelem_elem[lel];
1554 int nd = el_dof->RowSize(gel);
1555 int *gdofs = el_dof->GetRow(gel);
1557 for (
int j = 0; j < nd; j++)
1559 weights(gdofs[j]) = lext->
weights(ldofs[j]);
1572 for (
int i = 0; i < num_pieces; i++)
1579 for (
int lel = 0; lel < lext->
GetNE(); lel++)
1595 for (
int p = 0; p < GetNP(); p++)
1597 patchTopo->GetElementEdges(p, edges, oedge);
1599 for (
int i = 0; i < edges.
Size(); i++)
1601 edges[i] = edge_to_knot[edges[i]];
1604 edges[i] = -1 - edges[i];
1608 if ((Dimension() == 2 &&
1609 (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
1611 (Dimension() == 3 &&
1612 (edges[0] != edges[2] || edges[0] != edges[4] ||
1613 edges[0] != edges[6] || edges[1] != edges[3] ||
1614 edges[1] != edges[5] || edges[1] != edges[7] ||
1615 edges[8] != edges[9] || edges[8] != edges[10] ||
1616 edges[8] != edges[11])))
1618 cerr <<
"NURBSExtension::CheckPatch (patch = " << p
1619 <<
")\n Inconsistent edge-to-knot mapping!\n";
1623 if ((Dimension() == 2 &&
1624 (edges[0] < 0 || edges[1] < 0)) ||
1626 (Dimension() == 3 &&
1627 (edges[0] < 0 || edges[3] < 0 || edges[8] < 0)))
1629 cerr <<
"NURBSExtension::CheckPatch (patch = " << p
1630 <<
") : Bad orientation!\n";
1641 for (
int p = 0; p < GetNBP(); p++)
1643 patchTopo->GetBdrElementEdges(p, edges, oedge);
1645 for (
int i = 0; i < edges.
Size(); i++)
1647 edges[i] = edge_to_knot[edges[i]];
1650 edges[i] = -1 - edges[i];
1654 if ((Dimension() == 2 && (edges[0] < 0)) ||
1655 (Dimension() == 3 && (edges[0] < 0 || edges[1] < 0)))
1657 cerr <<
"NURBSExtension::CheckBdrPatch (boundary patch = "
1658 << p <<
") : Bad orientation!\n";
1670 patchTopo->GetElementEdges(p, edges, orient);
1672 if (Dimension() == 2)
1674 kv[0] = KnotVec(edges[0]);
1675 kv[1] = KnotVec(edges[1]);
1679 kv[0] = KnotVec(edges[0]);
1680 kv[1] = KnotVec(edges[3]);
1681 kv[2] = KnotVec(edges[8]);
1691 patchTopo->GetBdrElementEdges(p, edges, orient);
1693 if (Dimension() == 2)
1695 kv[0] = KnotVec(edges[0]);
1699 kv[0] = KnotVec(edges[0]);
1700 kv[1] = KnotVec(edges[1]);
1706 int nv = patchTopo->GetNV();
1707 int ne = patchTopo->GetNEdges();
1708 int nf = patchTopo->GetNFaces();
1709 int np = patchTopo->GetNE();
1710 int meshCounter, spaceCounter,
dim = Dimension();
1716 e_meshOffsets.SetSize(ne);
1717 f_meshOffsets.SetSize(nf);
1718 p_meshOffsets.SetSize(np);
1720 v_spaceOffsets.SetSize(nv);
1721 e_spaceOffsets.SetSize(ne);
1722 f_spaceOffsets.SetSize(nf);
1723 p_spaceOffsets.SetSize(np);
1726 for (meshCounter = 0; meshCounter < nv; meshCounter++)
1728 v_meshOffsets[meshCounter] = meshCounter;
1729 v_spaceOffsets[meshCounter] = meshCounter;
1731 spaceCounter = meshCounter;
1734 for (
int e = 0; e < ne; e++)
1736 e_meshOffsets[e] = meshCounter;
1737 e_spaceOffsets[e] = spaceCounter;
1738 meshCounter += KnotVec(e)->GetNE() - 1;
1739 spaceCounter += KnotVec(e)->GetNCP() - 2;
1743 for (
int f = 0; f < nf; f++)
1745 f_meshOffsets[f] = meshCounter;
1746 f_spaceOffsets[f] = spaceCounter;
1748 patchTopo->GetFaceEdges(f, edges, orient);
1751 (KnotVec(edges[0])->GetNE() - 1) *
1752 (KnotVec(edges[1])->GetNE() - 1);
1754 (KnotVec(edges[0])->GetNCP() - 2) *
1755 (KnotVec(edges[1])->GetNCP() - 2);
1759 for (
int p = 0; p < np; p++)
1761 p_meshOffsets[p] = meshCounter;
1762 p_spaceOffsets[p] = spaceCounter;
1764 patchTopo->GetElementEdges(p, edges, orient);
1769 (KnotVec(edges[0])->GetNE() - 1) *
1770 (KnotVec(edges[1])->GetNE() - 1);
1772 (KnotVec(edges[0])->GetNCP() - 2) *
1773 (KnotVec(edges[1])->GetNCP() - 2);
1778 (KnotVec(edges[0])->GetNE() - 1) *
1779 (KnotVec(edges[3])->GetNE() - 1) *
1780 (KnotVec(edges[8])->GetNE() - 1);
1782 (KnotVec(edges[0])->GetNCP() - 2) *
1783 (KnotVec(edges[3])->GetNCP() - 2) *
1784 (KnotVec(edges[8])->GetNCP() - 2);
1787 NumOfVertices = meshCounter;
1788 NumOfDofs = spaceCounter;
1793 int dim = Dimension();
1797 for (
int p = 0; p < GetNP(); p++)
1799 GetPatchKnotVectors(p, kv);
1801 int ne = kv[0]->GetNE();
1802 for (
int d = 1; d <
dim; d++)
1804 ne *= kv[d]->GetNE();
1807 NumOfElements += ne;
1813 int dim = Dimension() - 1;
1816 NumOfBdrElements = 0;
1817 for (
int p = 0; p < GetNBP(); p++)
1819 GetBdrPatchKnotVectors(p, kv);
1821 int ne = kv[0]->GetNE();
1822 for (
int d = 1; d <
dim; d++)
1824 ne *= kv[d]->GetNE();
1827 NumOfBdrElements += ne;
1835 if (Dimension() == 2)
1837 Get2DElementTopo(elements);
1841 Get3DElementTopo(elements);
1853 for (
int p = 0; p < GetNP(); p++)
1859 int patch_attr = patchTopo->GetAttribute(p);
1861 for (
int j = 0; j < ny; j++)
1863 for (
int i = 0; i < nx; i++)
1867 ind[0] = activeVert[p2g(i, j )];
1868 ind[1] = activeVert[p2g(i+1,j )];
1869 ind[2] = activeVert[p2g(i+1,j+1)];
1870 ind[3] = activeVert[p2g(i, j+1)];
1889 for (
int p = 0; p < GetNP(); p++)
1896 int patch_attr = patchTopo->GetAttribute(p);
1898 for (
int k = 0; k < nz; k++)
1900 for (
int j = 0; j < ny; j++)
1902 for (
int i = 0; i < nx; i++)
1906 ind[0] = activeVert[p2g(i, j, k)];
1907 ind[1] = activeVert[p2g(i+1,j, k)];
1908 ind[2] = activeVert[p2g(i+1,j+1,k)];
1909 ind[3] = activeVert[p2g(i, j+1,k)];
1911 ind[4] = activeVert[p2g(i, j, k+1)];
1912 ind[5] = activeVert[p2g(i+1,j, k+1)];
1913 ind[6] = activeVert[p2g(i+1,j+1,k+1)];
1914 ind[7] = activeVert[p2g(i, j+1,k+1)];
1916 elements[el] =
new Hexahedron(ind, patch_attr);
1930 if (Dimension() == 2)
1932 Get2DBdrElementTopo(boundary);
1936 Get3DBdrElementTopo(boundary);
1948 for (
int b = 0; b < GetNBP(); b++)
1953 int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
1955 for (
int i = 0; i < nx; i++)
1957 if (activeBdrElem[g_be])
1959 int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
1960 ind[0] = activeVert[p2g[_i ]];
1961 ind[1] = activeVert[p2g[_i+1]];
1963 boundary[l_be] =
new Segment(ind, bdr_patch_attr);
1979 for (
int b = 0; b < GetNBP(); b++)
1985 int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
1987 for (
int j = 0; j < ny; j++)
1989 int _j = (okv[1] >= 0) ? j : (ny - 1 - j);
1990 for (
int i = 0; i < nx; i++)
1992 if (activeBdrElem[g_be])
1994 int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
1995 ind[0] = activeVert[p2g(_i, _j )];
1996 ind[1] = activeVert[p2g(_i+1,_j )];
1997 ind[2] = activeVert[p2g(_i+1,_j+1)];
1998 ind[3] = activeVert[p2g(_i, _j+1)];
2011 activeDof.SetSize(GetNTotalDof());
2014 if (Dimension() == 2)
2016 Generate2DElementDofTable();
2020 Generate3DElementDofTable();
2023 NumOfActiveDofs = 0;
2024 for (
int d = 0; d < GetNTotalDof(); d++)
2028 activeDof[d] = NumOfActiveDofs;
2031 int *dof = el_dof->GetJ();
2032 int ndof = el_dof->Size_of_connections();
2033 for (
int i = 0; i < ndof; i++)
2035 dof[i] = activeDof[dof[i]] - 1;
2046 el_dof =
new Table(NumOfActiveElems, (Order + 1)*(Order + 1));
2047 el_to_patch.SetSize(NumOfActiveElems);
2048 el_to_IJK.SetSize(NumOfActiveElems, 2);
2050 for (
int p = 0; p < GetNP(); p++)
2055 for (
int j = 0; j < kv[1]->
GetNKS(); j++)
2057 if (kv[1]->isElement(j))
2059 for (
int i = 0; i < kv[0]->
GetNKS(); i++)
2061 if (kv[0]->isElement(i))
2065 int *dofs = el_dof->GetRow(el);
2067 for (
int jj = 0; jj <= Order; jj++)
2069 for (
int ii = 0; ii <= Order; ii++)
2071 dofs[idx] = p2g(i+ii,j+jj);
2072 activeDof[dofs[idx]] = 1;
2076 el_to_patch[el] = p;
2077 el_to_IJK(el,0) = i;
2078 el_to_IJK(el,1) = j;
2098 el_dof =
new Table(NumOfActiveElems, (Order + 1)*(Order + 1)*(Order + 1));
2099 el_to_patch.SetSize(NumOfActiveElems);
2100 el_to_IJK.SetSize(NumOfActiveElems, 3);
2102 for (
int p = 0; p < GetNP(); p++)
2107 for (
int k = 0; k < kv[2]->
GetNKS(); k++)
2109 if (kv[2]->isElement(k))
2111 for (
int j = 0; j < kv[1]->
GetNKS(); j++)
2113 if (kv[1]->isElement(j))
2115 for (
int i = 0; i < kv[0]->
GetNKS(); i++)
2117 if (kv[0]->isElement(i))
2122 int *dofs = el_dof->GetRow(el);
2123 for (
int kk = 0; kk <= Order; kk++)
2125 for (
int jj = 0; jj <= Order; jj++)
2127 for (
int ii = 0; ii <= Order; ii++)
2129 dofs[idx] = p2g(i+ii,j+jj,k+kk);
2130 activeDof[dofs[idx]] = 1;
2136 el_to_patch[el] = p;
2137 el_to_IJK(el,0) = i;
2138 el_to_IJK(el,1) = j;
2139 el_to_IJK(el,2) = k;
2155 if (Dimension() == 2)
2157 Generate2DBdrElementDofTable();
2161 Generate3DBdrElementDofTable();
2164 int *dof = bel_dof->GetJ();
2165 int ndof = bel_dof->Size_of_connections();
2166 for (
int i = 0; i < ndof; i++)
2168 dof[i] = activeDof[dof[i]] - 1;
2175 int lbe = 0, okv[1];
2179 bel_dof =
new Table(NumOfActiveBdrElems, Order + 1);
2180 bel_to_patch.SetSize(NumOfActiveBdrElems);
2181 bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
2183 for (
int b = 0; b < GetNBP(); b++)
2189 int nks0 = kv[0]->
GetNKS();
2190 for (
int i = 0; i < nks0; i++)
2192 if (kv[0]->isElement(i))
2194 if (activeBdrElem[gbe])
2196 int *dofs = bel_dof->GetRow(lbe);
2198 for (
int ii = 0; ii <= Order; ii++)
2200 dofs[idx] = p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)];
2203 bel_to_patch[lbe] = b;
2204 bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2216 int lbe = 0, okv[2];
2220 bel_dof =
new Table(NumOfActiveBdrElems, (Order + 1)*(Order + 1));
2221 bel_to_patch.SetSize(NumOfActiveBdrElems);
2222 bel_to_IJK.SetSize(NumOfActiveBdrElems, 2);
2224 for (
int b = 0; b < GetNBP(); b++)
2231 int nks0 = kv[0]->
GetNKS();
2232 int nks1 = kv[1]->
GetNKS();
2233 for (
int j = 0; j < nks1; j++)
2235 if (kv[1]->isElement(j))
2237 for (
int i = 0; i < nks0; i++)
2239 if (kv[0]->isElement(i))
2241 if (activeBdrElem[gbe])
2243 int *dofs = bel_dof->GetRow(lbe);
2245 for (
int jj = 0; jj <= Order; jj++)
2247 int _jj = (okv[1] >= 0) ? (j+jj) : (ny-j-jj);
2248 for (
int ii = 0; ii <= Order; ii++)
2250 int _ii = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
2251 dofs[idx] = p2g(_ii,_jj);
2255 bel_to_patch[lbe] = b;
2256 bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2257 bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
2271 for (
int gv = 0; gv < GetGNV(); gv++)
2272 if (activeVert[gv] >= 0)
2274 lvert_vert[activeVert[gv]] = gv;
2281 for (
int le = 0, ge = 0; ge < GetGNE(); ge++)
2284 lelem_elem[le++] = ge;
2296 NURBSFE->
SetIJK(el_to_IJK.GetRow(i));
2297 if (el_to_patch[i] != NURBSFE->
GetPatch())
2299 GetPatchKnotVectors(el_to_patch[i], NURBSFE->
KnotVectors());
2302 el_dof->GetRow(i, dofs);
2303 weights.GetSubVector(dofs, NURBSFE->
Weights());
2316 NURBSFE->
SetIJK(bel_to_IJK.GetRow(i));
2317 if (bel_to_patch[i] != NURBSFE->
GetPatch())
2319 GetBdrPatchKnotVectors(bel_to_patch[i], NURBSFE->
KnotVectors());
2320 NURBSFE->
SetPatch(bel_to_patch[i]);
2322 bel_dof->GetRow(i, dofs);
2323 weights.GetSubVector(dofs, NURBSFE->
Weights());
2333 if (patches.Size() == 0)
2335 GetPatchNets(Nodes);
2341 if (patches.Size() == 0) {
return; }
2343 SetSolutionVector(Nodes);
2349 if (patches.Size() == 0)
2350 mfem_error(
"NURBSExtension::SetKnotsFromPatches :"
2351 " No patches available!");
2355 for (
int p = 0; p < patches.Size(); p++)
2357 GetPatchKnotVectors(p, kv);
2359 for (
int i = 0; i < kv.
Size(); i++)
2361 *kv[i] = *patches[p]->GetKV(i);
2365 Order = knotVectors[0]->GetOrder();
2366 for (
int i = 0; i < NumOfKnotVectors; i++)
2368 if (Order != knotVectors[i]->GetOrder())
2370 " Variable orders are not supported!");
2378 NumOfActiveElems = NumOfElements;
2379 activeElem.SetSize(NumOfElements);
2382 GenerateActiveVertices();
2383 GenerateElementDofTable();
2384 GenerateActiveBdrElems();
2385 GenerateBdrElementDofTable();
2390 for (
int p = 0; p < patches.Size(); p++)
2392 patches[p]->DegreeElevate(t);
2398 for (
int p = 0; p < patches.Size(); p++)
2400 patches[p]->UniformRefinement();
2411 for (
int p = 0; p < patches.Size(); p++)
2413 patchTopo->GetElementEdges(p, edges, orient);
2417 pkv[0] = kv[KnotInd(edges[0])];
2418 pkv[1] = kv[KnotInd(edges[1])];
2422 pkv[0] = kv[KnotInd(edges[0])];
2423 pkv[1] = kv[KnotInd(edges[3])];
2424 pkv[2] = kv[KnotInd(edges[8])];
2427 patches[p]->KnotInsert(pkv);
2433 if (Dimension() == 2)
2435 Get2DPatchNets(coords);
2439 Get3DPatchNets(coords);
2448 patches.SetSize(GetNP());
2449 for (
int p = 0; p < GetNP(); p++)
2455 for (
int j = 0; j < kv[1]->GetNCP(); j++)
2457 for (
int i = 0; i < kv[0]->GetNCP(); i++)
2460 for (
int d = 0; d < 2; d++)
2462 Patch(i,j,d) = coords(l*2 + d)*weights(l);
2464 Patch(i,j,2) = weights(l);
2475 patches.SetSize(GetNP());
2476 for (
int p = 0; p < GetNP(); p++)
2482 for (
int k = 0; k < kv[2]->GetNCP(); k++)
2484 for (
int j = 0; j < kv[1]->GetNCP(); j++)
2486 for (
int i = 0; i < kv[0]->GetNCP(); i++)
2489 for (
int d = 0; d < 3; d++)
2491 Patch(i,j,k,d) = coords(l*3 + d)*weights(l);
2493 Patch(i,j,k,3) = weights(l);
2502 if (Dimension() == 2)
2504 Set2DSolutionVector(coords);
2508 Set3DSolutionVector(coords);
2517 weights.SetSize(GetNDof());
2518 for (
int p = 0; p < GetNP(); p++)
2523 for (
int j = 0; j < kv[1]->GetNCP(); j++)
2525 for (
int i = 0; i < kv[0]->GetNCP(); i++)
2528 for (
int d = 0; d < 2; d++)
2530 coords(l*2 + d) = Patch(i,j,d)/Patch(i,j,2);
2532 weights(l) = Patch(i,j,2);
2544 weights.SetSize(GetNDof());
2545 for (
int p = 0; p < GetNP(); p++)
2550 for (
int k = 0; k < kv[2]->GetNCP(); k++)
2552 for (
int j = 0; j < kv[1]->GetNCP(); j++)
2554 for (
int i = 0; i < kv[0]->GetNCP(); i++)
2557 for (
int d = 0; d < 3; d++)
2559 coords(l*3 + d) = Patch(i,j,k,d)/Patch(i,j,k,3);
2561 weights(l) = Patch(i,j,k,3);
2578 mfem_error(
"ParNURBSExtension::ParNURBSExtension :\n"
2579 " all elements in the parent must be active!");
2584 mfem_error(
"ParNURBSExtension::ParNURBSExtension :\n"
2585 " parent does not own the patch topology!");
2605 partitioning =
new int[
GetGNE()];
2606 for (
int i = 0; i <
GetGNE(); i++)
2608 partitioning[i] = part[i];
2610 SetActive(partitioning, active_bel);
2618 BuildGroups(partitioning, *serial_elem_dof);
2622 for (
int gel = 0, lel = 0; gel <
GetGNE(); gel++)
2628 int *gdofs = serial_elem_dof->
GetRow(gel);
2629 for (
int i = 0; i < ndofs; i++)
2640 : gtopo(par_parent->gtopo.GetComm())
2692 partitioning = NULL;
2695 if (par_parent->partitioning == NULL)
2696 mfem_error(
"ParNURBSExtension::ParNURBSExtension :\n"
2697 " parent ParNURBSExtension has no partitioning!");
2700 Table *serial_elem_dof = GetGlobalElementDofTable();
2701 BuildGroups(par_parent->partitioning, *serial_elem_dof);
2702 delete serial_elem_dof;
2705 Table *ParNURBSExtension::GetGlobalElementDofTable()
2709 return Get2DGlobalElementDofTable();
2713 return Get3DGlobalElementDofTable();
2717 Table *ParNURBSExtension::Get2DGlobalElementDofTable()
2725 for (
int p = 0; p <
GetNP(); p++)
2727 p2g.SetPatchDofMap(p, kv);
2730 for (
int j = 0; j < kv[1]->GetNKS(); j++)
2732 if (kv[1]->isElement(j))
2734 for (
int i = 0; i < kv[0]->GetNKS(); i++)
2736 if (kv[0]->isElement(i))
2738 int *dofs = gel_dof->GetRow(el);
2740 for (
int jj = 0; jj <=
Order; jj++)
2742 for (
int ii = 0; ii <=
Order; ii++)
2744 dofs[idx] = p2g(i+ii,j+jj);
2757 Table *ParNURBSExtension::Get3DGlobalElementDofTable()
2765 for (
int p = 0; p <
GetNP(); p++)
2767 p2g.SetPatchDofMap(p, kv);
2770 for (
int k = 0; k < kv[2]->GetNKS(); k++)
2772 if (kv[2]->isElement(k))
2774 for (
int j = 0; j < kv[1]->GetNKS(); j++)
2776 if (kv[1]->isElement(j))
2778 for (
int i = 0; i < kv[0]->GetNKS(); i++)
2780 if (kv[0]->isElement(i))
2782 int *dofs = gel_dof->GetRow(el);
2784 for (
int kk = 0; kk <=
Order; kk++)
2786 for (
int jj = 0; jj <=
Order; jj++)
2788 for (
int ii = 0; ii <=
Order; ii++)
2790 dofs[idx] = p2g(i+ii,j+jj,k+kk);
2806 void ParNURBSExtension::SetActive(
int *_partitioning,
2807 const Array<bool> &active_bel)
2813 for (
int i = 0; i <
GetGNE(); i++)
2814 if (_partitioning[i] == MyRank)
2822 for (
int i = 0; i <
GetGNBE(); i++)
2829 void ParNURBSExtension::BuildGroups(
int *_partitioning,
const Table &elem_dof)
2833 ListOfIntegerSets groups;
2838 for (
int i = 0; i < dof_proc.Size_of_connections(); i++)
2840 dof_proc.GetJ()[i] = _partitioning[dof_proc.GetJ()[i]];
2845 group.Recreate(1, &MyRank);
2846 groups.Insert(group);
2853 group.Recreate(dof_proc.RowSize(d), dof_proc.GetRow(d));
2864 void NURBSPatchMap::GetPatchKnotVectors(
int p, KnotVector *kv[])
2870 kv[0] = Ext->
KnotVec(edges[0]);
2871 kv[1] = Ext->
KnotVec(edges[1]);
2877 kv[0] = Ext->
KnotVec(edges[0]);
2878 kv[1] = Ext->
KnotVec(edges[3]);
2879 kv[2] = Ext->
KnotVec(edges[8]);
2884 void NURBSPatchMap::GetBdrPatchKnotVectors(
int p, KnotVector *kv[],
int *okv)
2888 kv[0] = Ext->
KnotVec(edges[0], oedge[0], &okv[0]);
2893 kv[1] = Ext->
KnotVec(edges[1], oedge[1], &okv[1]);
2903 GetPatchKnotVectors(p, kv);
2905 I = kv[0]->
GetNE() - 1;
2906 J = kv[1]->
GetNE() - 1;
2908 for (
int i = 0; i < verts.
Size(); i++)
2913 for (
int i = 0; i < edges.
Size(); i++)
2920 K = kv[2]->
GetNE() - 1;
2922 for (
int i = 0; i < faces.
Size(); i++)
2933 GetPatchKnotVectors(p, kv);
2938 for (
int i = 0; i < verts.
Size(); i++)
2943 for (
int i = 0; i < edges.
Size(); i++)
2952 for (
int i = 0; i < faces.
Size(); i++)
2963 GetBdrPatchKnotVectors(p, kv, okv);
2965 I = kv[0]->
GetNE() - 1;
2967 for (
int i = 0; i < verts.
Size(); i++)
2978 J = kv[1]->
GetNE() - 1;
2980 for (
int i = 0; i < edges.
Size(); i++)
2991 GetBdrPatchKnotVectors(p, kv, okv);
2995 for (
int i = 0; i < verts.
Size(); i++)
3008 for (
int i = 0; i < edges.
Size(); i++)
Array< KnotVector * > knotVectors
Abstract class for Finite Elements.
Array< int > f_meshOffsets
KnotVector * DegreeElevate(int t) const
void Get3DPatchNets(const Vector &Nodes)
void Set3DSolutionVector(Vector &Nodes)
void Create(ListOfIntegerSets &groups, int mpitag)
NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP)
Array2D< int > bel_to_IJK
int Size() const
Logical size of the array.
void SetCoordsFromPatches(Vector &Nodes)
KnotVector * KnotVec(int edge)
void swap(NURBSPatch *np)
void Get2DBdrElementTopo(Array< Element * > &boundary)
Class for grid function - Vector with associated FE space.
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
void DegreeElevate(int dir, int t)
NURBSPatch * Interpolate(NURBSPatch &p1, NURBSPatch &p2)
void SetPatchDofMap(int p, KnotVector *kv[])
void Generate2DElementDofTable()
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
void Generate3DElementDofTable()
void KnotInsert(int dir, const KnotVector &knot)
void MergeWeights(Mesh *mesh_array[], int num_pieces)
void SetSize(int s)
Resize the vector if the new size is different.
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
void GetElementFaces(int i, Array< int > &, Array< int > &) const
Return the indices and the orientations of all faces of element i.
void Print(std::ostream &out)
void Copy(Array ©) const
Create a copy of the current array.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
void Rotate3D(double normal[], double angle)
Array< int > e_meshOffsets
void GetElements()
Count the number of elements.
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
void GetPatchNets(const Vector &Nodes)
void Generate3DBdrElementDofTable()
Array< int > e_spaceOffsets
void GenerateBdrElementDofTable()
void GetElementVertices(int i, Array< int > &dofs) const
Returns the indices of the dofs of element i.
const KnotVector * GetKnotVector(int i) const
void Get2DElementTopo(Array< Element * > &elements)
void GetElementLocalToGlobal(Array< int > &lelem_elem)
Data type quadrilateral element.
Array< int > edge_to_knot
void GenerateActiveVertices()
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
void SetBdrPatchVertexMap(int p, KnotVector *kv[], int *okv)
void SwapDirections(int dir1, int dir2)
friend class NURBSPatchMap
void Get3DElementTopo(Array< Element * > &elements)
KnotVector()
Create KnotVector.
void Get3DBdrElementTopo(Array< Element * > &boundary)
static const int MaxOrder
Array< int > v_spaceOffsets
void UniformRefinement(Vector &newknots) const
void LoadFE(int i, const FiniteElement *FE)
void GetBdrElementTopo(Array< Element * > &boundary)
void GetElementTopo(Array< Element * > &elements)
void Generate2DBdrElementDofTable()
Data type hexahedron element.
void SetPatch(int p) const
void GetPatchKnotVectors(int p, Array< KnotVector * > &kv)
void GenerateActiveBdrElems()
Mesh * GetMesh() const
Returns the mesh.
void LoadPatchTopo(std::istream &input, Array< int > &edge_to_knot)
Read NURBS patch/macro-element mesh.
virtual ~NURBSExtension()
Destroy a NURBSExtension.
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
static void Get3DRotationMatrix(double n[], double angle, double r, DenseMatrix &T)
void SetKnotsFromPatches()
void CalcShape(Vector &shape, int i, double xi)
Array< int > f_spaceOffsets
Array< bool > activeBdrElem
FiniteElementSpace * FESpace()
void SetElement(int e) const
void GetBdrElementFace(int i, int *, int *) const
Return the index and the orientation of the face of bdr element i. (3D)
void Swap(Array< T > &, Array< T > &)
Array< int > v_meshOffsets
void mfem_error(const char *msg)
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
void GetBdrPatchKnotVectors(int p, Array< KnotVector * > &kv)
void SetSubVector(const Array< int > &dofs, const Vector &elemvect)
Array< int > p_meshOffsets
KnotVector & operator=(const KnotVector &kv)
void Set2DSolutionVector(Vector &Nodes)
Table * GetElementDofTable()
NURBSExtension * NURBSext
ParNURBSExtension(MPI_Comm comm, NURBSExtension *parent, int *partitioning, const Array< bool > &active_bel)
Array< KnotVector * > & KnotVectors() const
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
NURBSPatch * Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
Array< int > p_spaceOffsets
int findKnotSpan(double u) const
void GenerateElementDofTable()
static void skip_comment_lines(std::istream &is, const char comment_char)
void GetBdrElementVertices(int i, Array< int > &dofs) const
Returns the indices of the dofs of boundary element i.
void SetIJK(int *IJK) const
void SetBdrPatchDofMap(int p, KnotVector *kv[], int *okv)
void Difference(const KnotVector &kv, Vector &diff) const
void CalcDShape(Vector &grad, int i, double xi)
void DegreeElevate(int t)
void Get2DPatchNets(const Vector &Nodes)
void Mult(const double *x, double *y) const
Matrix vector multiplication.
void LoadBE(int i, const FiniteElement *BE)
void Print(std::ostream &out) const
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void PrintCharacteristics(std::ostream &out)
void SetPatchVertexMap(int p, KnotVector *kv[])
void FlipDirection(int dir)
void SetSolutionVector(Vector &Nodes)
Data type line segment element.
Array< int > bel_to_patch
void ConvertToPatches(const Vector &Nodes)
void KnotInsert(Array< KnotVector * > &kv)
void Print(std::ostream &out) const
int SetLoopDirection(int dir)