MFEM v2.0
|
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at 00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights 00003 // reserved. See file COPYRIGHT for details. 00004 // 00005 // This file is part of the MFEM library. For more information and source code 00006 // availability see http://mfem.googlecode.com. 00007 // 00008 // MFEM is free software; you can redistribute it and/or modify it under the 00009 // terms of the GNU Lesser General Public License (as published by the Free 00010 // Software Foundation) version 2.1 dated February 1999. 00011 00012 #include "../fem/fem.hpp" 00013 00014 00015 const int KnotVector::MaxOrder = 10; 00016 00017 KnotVector::KnotVector(istream &input) 00018 { 00019 input >> Order >> NumOfControlPoints; 00020 knot.Load(input, NumOfControlPoints + Order + 1); 00021 GetElements(); 00022 } 00023 00024 KnotVector::KnotVector(int Order_, int NCP) 00025 { 00026 Order = Order_; 00027 NumOfControlPoints = NCP; 00028 knot.SetSize(NumOfControlPoints + Order + 1); 00029 00030 knot = -1.; 00031 } 00032 00033 KnotVector &KnotVector::operator=(const KnotVector &kv) 00034 { 00035 Order = kv.Order; 00036 NumOfControlPoints = kv.NumOfControlPoints; 00037 NumOfElements = kv.NumOfElements; 00038 knot = kv.knot; 00039 // alternatively, re-compute NumOfElements 00040 // GetElements(); 00041 00042 return *this; 00043 } 00044 00045 KnotVector *KnotVector::DegreeElevate(int t) const 00046 { 00047 if (t < 0) 00048 mfem_error("KnotVector::DegreeElevate :\n" 00049 " Parent KnotVector order higher than child"); 00050 00051 int nOrder = Order + t; 00052 KnotVector *newkv = new KnotVector(nOrder, GetNCP() + t); 00053 00054 for (int i = 0; i <= nOrder; i++) 00055 { 00056 (*newkv)[i] = knot(0); 00057 } 00058 for (int i = nOrder + 1; i < newkv->GetNCP(); i++) 00059 { 00060 (*newkv)[i] = knot(i - t); 00061 } 00062 for (int i = 0; i <= nOrder; i++) 00063 { 00064 (*newkv)[newkv->GetNCP() + i] = knot(knot.Size()-1); 00065 } 00066 00067 newkv->GetElements(); 00068 00069 return newkv; 00070 } 00071 00072 void KnotVector::UniformRefinement(Vector &newknots) const 00073 { 00074 newknots.SetSize(NumOfElements); 00075 int j = 0; 00076 for (int i = 0; i < knot.Size()-1; i++) 00077 { 00078 if (knot(i) != knot(i+1)) 00079 { 00080 newknots(j) = 0.5*(knot(i) + knot(i+1)); 00081 j++; 00082 } 00083 } 00084 } 00085 00086 void KnotVector::GetElements() 00087 { 00088 NumOfElements = 0; 00089 for (int i = Order; i < NumOfControlPoints; i++) 00090 { 00091 if (knot(i) != knot(i+1)) 00092 { 00093 NumOfElements++; 00094 } 00095 } 00096 } 00097 00098 void KnotVector::Flip() 00099 { 00100 double apb = knot(0) + knot(knot.Size()-1); 00101 00102 int ns = (NumOfControlPoints - Order)/2; 00103 for (int i = 1; i <= ns; i++) 00104 { 00105 double tmp = apb - knot(Order + i); 00106 knot(Order + i) = apb - knot(NumOfControlPoints - i); 00107 knot(NumOfControlPoints - i) = tmp; 00108 } 00109 } 00110 00111 void KnotVector::Print(ostream &out) const 00112 { 00113 out << Order << ' ' << NumOfControlPoints << ' '; 00114 knot.Print(out, knot.Size()); 00115 } 00116 00117 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller 00118 void KnotVector::CalcShape(Vector &shape, int i, double xi) 00119 { 00120 int p = Order; 00121 int ip = (i >= 0) ? (i + p) : (-1 - i + p); 00122 double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), saved, tmp; 00123 double left[MaxOrder+1], right[MaxOrder+1]; 00124 00125 #ifdef MFEM_DEBUG 00126 if (p > MaxOrder) 00127 mfem_error("KnotVector::CalcShape : Order > MaxOrder!"); 00128 #endif 00129 00130 shape(0) = 1.; 00131 for (int j = 1; j <= p; ++j) 00132 { 00133 left[j] = u - knot(ip+1-j); 00134 right[j] = knot(ip+j) - u; 00135 saved = 0.; 00136 for (int r = 0; r < j; ++r) 00137 { 00138 tmp = shape(r)/(right[r+1] + left[j-r]); 00139 shape(r) = saved + right[r+1]*tmp; 00140 saved = left[j-r]*tmp; 00141 } 00142 shape(j) = saved; 00143 } 00144 } 00145 00146 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller 00147 void KnotVector::CalcDShape(Vector &grad, int i, double xi) 00148 { 00149 int p = Order, rk, pk; 00150 int ip = (i >= 0) ? (i + p) : (-1 - i + p); 00151 double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), temp, saved, d; 00152 double ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1], right[MaxOrder+1]; 00153 00154 #ifdef MFEM_DEBUG 00155 if (p > MaxOrder) 00156 mfem_error("KnotVector::CalcDShape : Order > MaxOrder!"); 00157 #endif 00158 00159 ndu[0][0] = 1.0; 00160 for (int j = 1; j <= p; j++) 00161 { 00162 left[j] = u - knot(ip-j+1); 00163 right[j] = knot(ip+j) - u; 00164 saved = 0.0; 00165 for (int r = 0; r < j; r++) 00166 { 00167 ndu[j][r] = right[r+1] + left[j-r]; 00168 temp = ndu[r][j-1]/ndu[j][r]; 00169 ndu[r][j] = saved + right[r+1]*temp; 00170 saved = left[j-r]*temp; 00171 } 00172 ndu[j][j] = saved; 00173 } 00174 00175 for (int r = 0; r <= p; ++r) 00176 { 00177 d = 0.0; 00178 rk = r-1; 00179 pk = p-1; 00180 if (r >= 1) 00181 { 00182 d = ndu[rk][pk]/ndu[p][rk]; 00183 } 00184 if (r <= pk) 00185 { 00186 d -= ndu[r][pk]/ndu[p][r]; 00187 } 00188 grad(r) = d; 00189 } 00190 00191 if (i >= 0) 00192 grad *= p*(knot(ip+1) - knot(ip)); 00193 else 00194 grad *= p*(knot(ip) - knot(ip+1)); 00195 } 00196 00197 int KnotVector::findKnotSpan(double u) const 00198 { 00199 int low, mid, high; 00200 00201 if (u == knot(NumOfControlPoints+Order)) 00202 { 00203 mid = NumOfControlPoints; 00204 } 00205 else 00206 { 00207 low = Order; 00208 high = NumOfControlPoints + 1; 00209 mid = (low + high)/2; 00210 while ( (u < knot(mid-1)) || (u > knot(mid)) ) 00211 { 00212 if (u < knot(mid-1)) 00213 { 00214 high = mid; 00215 } 00216 else 00217 { 00218 low = mid; 00219 } 00220 mid = (low + high)/2; 00221 } 00222 } 00223 return mid; 00224 } 00225 00226 void KnotVector::Difference(const KnotVector &kv, Vector &diff) const 00227 { 00228 if (Order != kv.GetOrder()) 00229 { 00230 mfem_error("KnotVector::Difference :\n" 00231 " Can not compare knot vectors with different orders!"); 00232 } 00233 00234 int s = kv.Size() - Size(); 00235 if (s < 0) 00236 { 00237 kv.Difference(*this, diff); 00238 return; 00239 } 00240 00241 diff.SetSize(s); 00242 00243 s = 0; 00244 int i = 0; 00245 for (int j = 0; j < kv.Size(); j++) 00246 { 00247 if (knot(i) == kv[j]) 00248 { 00249 i++; 00250 } 00251 else 00252 { 00253 diff(s) = kv[j]; 00254 s++; 00255 } 00256 } 00257 } 00258 00259 00260 void NURBSPatch::init(int dim_) 00261 { 00262 Dim = dim_; 00263 sd = nd = -1; 00264 00265 if (kv.Size() == 2) 00266 { 00267 ni = kv[0]->GetNCP(); 00268 nj = kv[1]->GetNCP(); 00269 nk = -1; 00270 00271 data = new double[ni*nj*Dim]; 00272 00273 #ifdef MFEM_DEBUG 00274 for (int i = 0; i < ni*nj*Dim; i++) 00275 data[i] = -999.99; 00276 #endif 00277 } 00278 else if (kv.Size() == 3) 00279 { 00280 ni = kv[0]->GetNCP(); 00281 nj = kv[1]->GetNCP(); 00282 nk = kv[2]->GetNCP(); 00283 00284 data = new double[ni*nj*nk*Dim]; 00285 00286 #ifdef MFEM_DEBUG 00287 for (int i = 0; i < ni*nj*nk*Dim; i++) 00288 data[i] = -999.99; 00289 #endif 00290 } 00291 else 00292 { 00293 mfem_error("NURBSPatch::init : Wrond dimension of knotvectors!"); 00294 } 00295 } 00296 00297 NURBSPatch::NURBSPatch(istream &input) 00298 { 00299 int pdim, dim, size = 1; 00300 string ident; 00301 00302 input >> ws >> ident >> pdim; // knotvectors 00303 kv.SetSize(pdim); 00304 for (int i = 0; i < pdim; i++) 00305 { 00306 kv[i] = new KnotVector(input); 00307 size *= kv[i]->GetNCP(); 00308 } 00309 00310 input >> ws >> ident >> dim; // dimension 00311 init(dim + 1); 00312 00313 input >> ws >> ident; // controlpoints (homogeneous coordinates) 00314 for (int j = 0, i = 0; i < size; i++) 00315 for (int d = 0; d <= dim; d++, j++) 00316 input >> data[j]; 00317 } 00318 00319 NURBSPatch::NURBSPatch(KnotVector *kv0, KnotVector *kv1, int dim_) 00320 { 00321 kv.SetSize(2); 00322 kv[0] = new KnotVector(*kv0); 00323 kv[1] = new KnotVector(*kv1); 00324 init(dim_); 00325 } 00326 00327 NURBSPatch::NURBSPatch(KnotVector *kv0, KnotVector *kv1, KnotVector *kv2, 00328 int dim_) 00329 { 00330 kv.SetSize(3); 00331 kv[0] = new KnotVector(*kv0); 00332 kv[1] = new KnotVector(*kv1); 00333 kv[2] = new KnotVector(*kv2); 00334 init(dim_); 00335 } 00336 00337 NURBSPatch::NURBSPatch(Array<KnotVector *> &kv_, int dim_) 00338 { 00339 kv.SetSize(kv_.Size()); 00340 for (int i = 0; i < kv.Size(); i++) 00341 kv[i] = new KnotVector(*kv_[i]); 00342 init(dim_); 00343 } 00344 00345 NURBSPatch::NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP) 00346 { 00347 kv.SetSize(parent->kv.Size()); 00348 for (int i = 0; i < kv.Size(); i++) 00349 if (i != dir) 00350 kv[i] = new KnotVector(*parent->kv[i]); 00351 else 00352 kv[i] = new KnotVector(Order, NCP); 00353 init(parent->Dim); 00354 } 00355 00356 void NURBSPatch::swap(NURBSPatch *np) 00357 { 00358 if (data != NULL) 00359 delete [] data; 00360 00361 for (int i = 0; i < kv.Size(); i++) 00362 { 00363 if (kv[i]) delete kv[i]; 00364 } 00365 00366 data = np->data; 00367 np->kv.Copy(kv); 00368 00369 ni = np->ni; 00370 nj = np->nj; 00371 nk = np->nk; 00372 Dim = np->Dim; 00373 00374 np->data = NULL; 00375 np->kv.SetSize(0); 00376 00377 delete np; 00378 } 00379 00380 NURBSPatch::~NURBSPatch() 00381 { 00382 if (data != NULL) 00383 delete [] data; 00384 00385 for (int i = 0; i < kv.Size(); i++) 00386 { 00387 if (kv[i]) delete kv[i]; 00388 } 00389 } 00390 00391 void NURBSPatch::Print(ostream &out) 00392 { 00393 int size = 1; 00394 00395 out << "knotvectors\n" << kv.Size() << '\n'; 00396 for (int i = 0; i < kv.Size(); i++) 00397 { 00398 kv[i]->Print(out); 00399 size *= kv[i]->GetNCP(); 00400 } 00401 00402 out << "\ndimension\n" << Dim - 1 00403 << "\n\ncontrolpoints\n"; 00404 for (int j = 0, i = 0; i < size; i++) 00405 { 00406 out << data[j++]; 00407 for (int d = 1; d < Dim; d++) 00408 out << ' ' << data[j++]; 00409 out << '\n'; 00410 } 00411 } 00412 00413 int NURBSPatch::SetLoopDirection(int dir) 00414 { 00415 if (nk == -1) 00416 { 00417 if (dir == 0) 00418 { 00419 sd = Dim; 00420 nd = ni; 00421 00422 return nj*Dim; 00423 } 00424 else if (dir == 1) 00425 { 00426 sd = ni*Dim; 00427 nd = nj; 00428 00429 return ni*Dim; 00430 } 00431 else 00432 { 00433 cerr << "NURBSPatch::SetLoopDirection :\n" 00434 " Direction error in 2D patch, dir = " << dir << '\n'; 00435 mfem_error(); 00436 } 00437 } 00438 else 00439 { 00440 if (dir == 0) 00441 { 00442 sd = Dim; 00443 nd = ni; 00444 00445 return nj*nk*Dim; 00446 } 00447 else if (dir == 1) 00448 { 00449 sd = ni*Dim; 00450 nd = nj; 00451 00452 return ni*nk*Dim; 00453 } 00454 else if (dir == 2) 00455 { 00456 sd = ni*nj*Dim; 00457 nd = nk; 00458 00459 return ni*nj*Dim; 00460 } 00461 else 00462 { 00463 cerr << "NURBSPatch::SetLoopDirection :\n" 00464 " Direction error in 3D patch, dir = " << dir << '\n'; 00465 mfem_error(); 00466 } 00467 } 00468 00469 return -1; 00470 } 00471 00472 void NURBSPatch::UniformRefinement() 00473 { 00474 Vector newknots; 00475 for (int dir = 0; dir < kv.Size(); dir++) 00476 { 00477 kv[dir]->UniformRefinement(newknots); 00478 KnotInsert(dir, newknots); 00479 } 00480 } 00481 00482 void NURBSPatch::KnotInsert(Array<KnotVector *> &newkv) 00483 { 00484 for (int dir = 0; dir < kv.Size(); dir++) 00485 { 00486 KnotInsert(dir, *newkv[dir]); 00487 } 00488 } 00489 00490 void NURBSPatch::KnotInsert(int dir, const KnotVector &newkv) 00491 { 00492 if (dir >= kv.Size() || dir < 0) 00493 mfem_error("NURBSPatch::KnotInsert : Incorrect direction!"); 00494 00495 int t = newkv.GetOrder() - kv[dir]->GetOrder(); 00496 00497 if (t > 0) 00498 { 00499 DegreeElevate(dir, t); 00500 } 00501 else if (t < 0) 00502 { 00503 mfem_error("NURBSPatch::KnotInsert : Incorrect order!"); 00504 } 00505 00506 Vector diff; 00507 GetKV(dir)->Difference(newkv, diff); 00508 if (diff.Size() > 0) 00509 { 00510 KnotInsert(dir, diff); 00511 } 00512 } 00513 00514 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller 00515 void NURBSPatch::KnotInsert(int dir, const Vector &knot) 00516 { 00517 if (dir >= kv.Size() || dir < 0) 00518 mfem_error("NURBSPatch::KnotInsert : Incorrect direction!"); 00519 00520 NURBSPatch &oldp = *this; 00521 KnotVector &oldkv = *kv[dir]; 00522 00523 NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder(), 00524 oldkv.GetNCP() + knot.Size()); 00525 NURBSPatch &newp = *newpatch; 00526 KnotVector &newkv = *newp.GetKV(dir); 00527 00528 int size = oldp.SetLoopDirection(dir); 00529 if (size != newp.SetLoopDirection(dir)) 00530 mfem_error("NURBSPatch::KnotInsert : Size mismatch!"); 00531 00532 int rr = knot.Size() - 1; 00533 int a = oldkv.findKnotSpan(knot(0)) - 1; 00534 int b = oldkv.findKnotSpan(knot(rr)) - 1; 00535 int pl = oldkv.GetOrder(); 00536 int ml = oldkv.GetNCP(); 00537 00538 for (int j = 0; j <= a; j++) 00539 { 00540 newkv[j] = oldkv[j]; 00541 } 00542 for (int j = b+pl; j <= ml+pl; j++) 00543 { 00544 newkv[j+rr+1] = oldkv[j]; 00545 } 00546 for (int k = 0; k <= (a-pl); k++) 00547 { 00548 for (int ll = 0; ll < size; ll++) 00549 newp(k,ll) = oldp(k,ll); 00550 } 00551 for (int k = (b-1); k < ml; k++) 00552 { 00553 for (int ll = 0; ll < size; ll++) 00554 newp(k+rr+1,ll) = oldp(k,ll); 00555 } 00556 00557 int i = b+pl-1; 00558 int k = b+pl+rr; 00559 00560 for (int j = rr; j >= 0; j--) 00561 { 00562 while ( (knot(j) <= oldkv[i]) && (i > a) ) 00563 { 00564 newkv[k] = oldkv[i]; 00565 for (int ll = 0; ll < size; ll++) 00566 newp(k-pl-1,ll) = oldp(i-pl-1,ll); 00567 00568 k--; 00569 i--; 00570 } 00571 00572 for (int ll = 0; ll < size; ll++) 00573 newp(k-pl-1,ll) = newp(k-pl,ll); 00574 00575 for (int l = 1; l <= pl; l++) 00576 { 00577 int ind = k-pl+l; 00578 double alfa = newkv[k+l] - knot(j); 00579 if (fabs(alfa) == 0.0) 00580 { 00581 for (int ll = 0; ll < size; ll++) 00582 newp(ind-1,ll) = newp(ind,ll); 00583 } 00584 else 00585 { 00586 alfa = alfa/(newkv[k+l] - oldkv[i-pl+l]); 00587 for (int ll = 0; ll < size; ll++) 00588 newp(ind-1,ll) = alfa*newp(ind-1,ll) + (1.0-alfa)*newp(ind,ll); 00589 } 00590 } 00591 00592 newkv[k] = knot(j); 00593 k--; 00594 } 00595 00596 newkv.GetElements(); 00597 00598 swap(newpatch); 00599 } 00600 00601 void NURBSPatch::DegreeElevate(int t) 00602 { 00603 for (int dir = 0; dir < kv.Size(); dir++) 00604 { 00605 DegreeElevate(dir, t); 00606 } 00607 } 00608 00609 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller 00610 void NURBSPatch::DegreeElevate(int dir, int t) 00611 { 00612 if (dir >= kv.Size() || dir < 0) 00613 mfem_error("NURBSPatch::DegreeElevate : Incorrect direction!"); 00614 00615 int i, j, k, kj, mpi, mul, mh, kind, cind, first, last; 00616 int r, a, b, oldr, save, s, tr, lbz, rbz, l; 00617 double inv, ua, ub, numer, alf, den, bet, gam; 00618 00619 NURBSPatch &oldp = *this; 00620 KnotVector &oldkv = *kv[dir]; 00621 00622 NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder() + t, 00623 oldkv.GetNCP() + oldkv.GetNE()*t); 00624 NURBSPatch &newp = *newpatch; 00625 KnotVector &newkv = *newp.GetKV(dir); 00626 00627 int size = oldp.SetLoopDirection(dir); 00628 if (size != newp.SetLoopDirection(dir)) 00629 mfem_error("NURBSPatch::DegreeElevate : Size mismatch!"); 00630 00631 int p = oldkv.GetOrder(); 00632 int n = oldkv.GetNCP()-1; 00633 00634 DenseMatrix bezalfs (p+t+1, p+1); 00635 DenseMatrix bpts (p+1, size); 00636 DenseMatrix ebpts (p+t+1, size); 00637 DenseMatrix nextbpts(p-1, size); 00638 Vector alphas (p-1); 00639 00640 int m = n + p + 1; 00641 int ph = p + t; 00642 int ph2 = ph/2; 00643 00644 { 00645 Array2D<int> binom(ph+1, ph+1); 00646 for (i = 0; i <= ph; i++) 00647 { 00648 binom(i,0) = binom(i,i) = 1; 00649 for (j = 1; j < i; j++) 00650 binom(i,j) = binom(i-1,j) + binom(i-1,j-1); 00651 } 00652 00653 bezalfs(0,0) = 1.0; 00654 bezalfs(ph,p) = 1.0; 00655 00656 for (i = 1; i <= ph2; i++) 00657 { 00658 inv = 1.0/binom(ph,i); 00659 mpi = min(p,i); 00660 for (j = max(0,i-t); j <= mpi; j++) 00661 { 00662 bezalfs(i,j) = inv*binom(p,j)*binom(t,i-j); 00663 } 00664 } 00665 } 00666 00667 for (i = ph2+1; i < ph; i++) 00668 { 00669 mpi = min(p,i); 00670 for (j = max(0,i-t); j <= mpi; j++) 00671 { 00672 bezalfs(i,j) = bezalfs(ph-i,p-j); 00673 } 00674 } 00675 00676 mh = ph; 00677 kind = ph + 1; 00678 r = -1; 00679 a = p; 00680 b = p + 1; 00681 cind = 1; 00682 ua = oldkv[0]; 00683 for (l = 0; l < size; l++) 00684 newp(0,l) = oldp(0,l); 00685 for (i = 0; i <= ph; i++) 00686 newkv[i] = ua; 00687 00688 for (i = 0; i <= p; i++) 00689 { 00690 for (l = 0; l < size; l++) 00691 bpts(i,l) = oldp(i,l); 00692 } 00693 00694 while (b < m) 00695 { 00696 i = b; 00697 while (b < m && oldkv[b] == oldkv[b+1]) b++; 00698 00699 mul = b-i+1; 00700 00701 mh = mh + mul + t; 00702 ub = oldkv[b]; 00703 oldr = r; 00704 r = p-mul; 00705 if (oldr > 0) lbz = (oldr+2)/2; 00706 else lbz = 1; 00707 00708 if (r > 0) rbz = ph-(r+1)/2; 00709 else rbz = ph; 00710 00711 if (r > 0) 00712 { 00713 numer = ub - ua; 00714 for (k = p ; k > mul; k--) 00715 alphas[k-mul-1] = numer/(oldkv[a+k]-ua); 00716 00717 for (j = 1; j <= r; j++) 00718 { 00719 save = r-j; 00720 s = mul+j; 00721 for (k = p; k >= s; k--) 00722 { 00723 for (l = 0; l < size; l++) 00724 bpts(k,l) = (alphas[k-s]*bpts(k,l) + 00725 (1.0-alphas[k-s])*bpts(k-1,l)); 00726 } 00727 for (l = 0; l < size; l++) 00728 nextbpts(save,l) = bpts(p,l); 00729 } 00730 } 00731 00732 for (i = lbz; i <= ph; i++) 00733 { 00734 for (l = 0; l < size; l++) 00735 ebpts(i,l) = 0.0; 00736 mpi = min(p,i); 00737 for (j = max(0,i-t); j <= mpi; j++) 00738 { 00739 for (l = 0; l < size; l++) 00740 ebpts(i,l) += bezalfs(i,j)*bpts(j,l); 00741 } 00742 } 00743 00744 if (oldr > 1) 00745 { 00746 first = kind-2; 00747 last = kind; 00748 den = ub-ua; 00749 bet = (ub-newkv[kind-1])/den; 00750 00751 for (tr = 1; tr < oldr; tr++) 00752 { 00753 i = first; 00754 j = last; 00755 kj = j-kind+1; 00756 while (j-i > tr) 00757 { 00758 if (i < cind) 00759 { 00760 alf = (ub-newkv[i])/(ua-newkv[i]); 00761 for (l = 0; l < size; l++) 00762 newp(i,l) = alf*newp(i,l)-(1.0-alf)*newp(i-1,l); 00763 } 00764 if (j >= lbz) 00765 { 00766 if ((j-tr) <= (kind-ph+oldr)) 00767 { 00768 gam = (ub-newkv[j-tr])/den; 00769 for (l = 0; l < size; l++) 00770 ebpts(kj,l) = gam*ebpts(kj,l) + (1.0-gam)*ebpts(kj+1,l); 00771 } 00772 else 00773 { 00774 for (l = 0; l < size; l++) 00775 ebpts(kj,l) = bet*ebpts(kj,l) + (1.0-bet)*ebpts(kj+1,l); 00776 } 00777 } 00778 i = i+1; 00779 j = j-1; 00780 kj = kj-1; 00781 } 00782 first--; 00783 last++; 00784 } 00785 } 00786 00787 if (a != p) 00788 { 00789 for (i = 0; i < (ph-oldr); i++) 00790 { 00791 newkv[kind] = ua; 00792 kind = kind+1; 00793 } 00794 } 00795 for (j = lbz; j <= rbz; j++) 00796 { 00797 for (l = 0; l < size; l++) 00798 newp(cind,l) = ebpts(j,l); 00799 cind = cind +1; 00800 } 00801 00802 if (b < m) 00803 { 00804 for (j = 0; j <r; j++) 00805 for (l = 0; l < size; l++) 00806 bpts(j,l) = nextbpts(j,l); 00807 00808 for (j = r; j <= p; j++) 00809 for (l = 0; l < size; l++) 00810 bpts(j,l) = oldp(b-p+j,l); 00811 00812 a = b; 00813 b = b+1; 00814 ua = ub; 00815 } 00816 else 00817 { 00818 for (i = 0; i <= ph; i++) 00819 newkv[kind+i] = ub; 00820 } 00821 } 00822 newkv.GetElements(); 00823 00824 swap(newpatch); 00825 } 00826 00827 void NURBSPatch::FlipDirection(int dir) 00828 { 00829 int size = SetLoopDirection(dir); 00830 00831 for (int id = 0; id < nd/2; id++) 00832 for (int i = 0; i < size; i++) 00833 Swap<double>((*this)(id,i), (*this)(nd-1-id,i)); 00834 kv[dir]->Flip(); 00835 } 00836 00837 void NURBSPatch::SwapDirections(int dir1, int dir2) 00838 { 00839 if (abs(dir1-dir2) == 2) 00840 mfem_error("NURBSPatch::SwapDirections :" 00841 " directions 0 and 2 are not supported!"); 00842 00843 Array<KnotVector *> nkv; 00844 00845 kv.Copy(nkv); 00846 Swap<KnotVector *>(nkv[dir1], nkv[dir2]); 00847 NURBSPatch *newpatch = new NURBSPatch(nkv, Dim); 00848 00849 int size = SetLoopDirection(dir1); 00850 newpatch->SetLoopDirection(dir2); 00851 00852 for (int id = 0; id < nd; id++) 00853 for (int i = 0; i < size; i++) 00854 (*newpatch)(id,i) = (*this)(id,i); 00855 00856 swap(newpatch); 00857 } 00858 00859 void NURBSPatch::Get3DRotationMatrix(double n[], double angle, double r, 00860 DenseMatrix &T) 00861 { 00862 double c, s, c1; 00863 double l2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2]; 00864 double l = sqrt(l2); 00865 00866 if (fabs(angle) == M_PI_2) 00867 { 00868 s = r*copysign(1., angle); 00869 c = 0.; 00870 c1 = -1.; 00871 } 00872 else if (fabs(angle) == M_PI) 00873 { 00874 s = 0.; 00875 c = -r; 00876 c1 = c - 1.; 00877 } 00878 else 00879 { 00880 s = r*sin(angle); 00881 c = r*cos(angle); 00882 c1 = c - 1.; 00883 } 00884 00885 T.SetSize(3); 00886 00887 T(0,0) = (n[0]*n[0] + (n[1]*n[1] + n[2]*n[2])*c)/l2; 00888 T(0,1) = -(n[0]*n[1]*c1)/l2 - (n[2]*s)/l; 00889 T(0,2) = -(n[0]*n[2]*c1)/l2 + (n[1]*s)/l; 00890 T(1,0) = -(n[0]*n[1]*c1)/l2 + (n[2]*s)/l; 00891 T(1,1) = (n[1]*n[1] + (n[0]*n[0] + n[2]*n[2])*c)/l2; 00892 T(1,2) = -(n[1]*n[2]*c1)/l2 - (n[0]*s)/l; 00893 T(2,0) = -(n[0]*n[2]*c1)/l2 - (n[1]*s)/l; 00894 T(2,1) = -(n[1]*n[2]*c1)/l2 + (n[0]*s)/l; 00895 T(2,2) = (n[2]*n[2] + (n[0]*n[0] + n[1]*n[1])*c)/l2; 00896 } 00897 00898 void NURBSPatch::Rotate3D(double n[], double angle) 00899 { 00900 if (Dim != 4) 00901 mfem_error("NURBSPatch::Rotate3D : not a NURBSPatch in 3D!"); 00902 00903 DenseMatrix T(3); 00904 Vector x(3), y(NULL, 3); 00905 00906 Get3DRotationMatrix(n, angle, 1., T); 00907 00908 int size = 1; 00909 for (int i = 0; i < kv.Size(); i++) 00910 size *= kv[i]->GetNCP(); 00911 00912 for (int i = 0; i < size; i++) 00913 { 00914 y.SetData(data + i*Dim); 00915 x = y; 00916 T.Mult(x, y); 00917 } 00918 } 00919 00920 int NURBSPatch::MakeUniformDegree() 00921 { 00922 int maxd = -1; 00923 00924 for (int dir = 0; dir < kv.Size(); dir++) 00925 if (maxd < kv[dir]->GetOrder()) 00926 maxd = kv[dir]->GetOrder(); 00927 00928 for (int dir = 0; dir < kv.Size(); dir++) 00929 if (maxd > kv[dir]->GetOrder()) 00930 DegreeElevate(dir, maxd - kv[dir]->GetOrder()); 00931 00932 return maxd; 00933 } 00934 00935 NURBSPatch *Interpolate(NURBSPatch &p1, NURBSPatch &p2) 00936 { 00937 if (p1.kv.Size() != p2.kv.Size() || p1.Dim != p2.Dim) 00938 mfem_error("Interpolate(NURBSPatch &, NURBSPatch &)"); 00939 00940 int size = 1, dim = p1.Dim; 00941 Array<KnotVector *> kv(p1.kv.Size() + 1); 00942 00943 for (int i = 0; i < p1.kv.Size(); i++) 00944 { 00945 if (p1.kv[i]->GetOrder() < p2.kv[i]->GetOrder()) 00946 { 00947 p1.KnotInsert(i, *p2.kv[i]); 00948 p2.KnotInsert(i, *p1.kv[i]); 00949 } 00950 else 00951 { 00952 p2.KnotInsert(i, *p1.kv[i]); 00953 p1.KnotInsert(i, *p2.kv[i]); 00954 } 00955 kv[i] = p1.kv[i]; 00956 size *= kv[i]->GetNCP(); 00957 } 00958 00959 kv.Last() = new KnotVector(1, 2); 00960 KnotVector &nkv = *kv.Last(); 00961 nkv[0] = nkv[1] = 0.0; 00962 nkv[2] = nkv[3] = 1.0; 00963 nkv.GetElements(); 00964 00965 NURBSPatch *patch = new NURBSPatch(kv, dim); 00966 delete kv.Last(); 00967 00968 for (int i = 0; i < size; i++) 00969 { 00970 for (int d = 0; d < dim; d++) 00971 { 00972 patch->data[i*dim+d] = p1.data[i*dim+d]; 00973 patch->data[(i+size)*dim+d] = p2.data[i*dim+d]; 00974 } 00975 } 00976 00977 return patch; 00978 } 00979 00980 NURBSPatch *Revolve3D(NURBSPatch &patch, double n[], double ang, int times) 00981 { 00982 if (patch.Dim != 4) 00983 mfem_error("Revolve3D(NURBSPatch &, double [], double)"); 00984 00985 int size = 1, ns; 00986 Array<KnotVector *> nkv(patch.kv.Size() + 1); 00987 00988 for (int i = 0; i < patch.kv.Size(); i++) 00989 { 00990 nkv[i] = patch.kv[i]; 00991 size *= nkv[i]->GetNCP(); 00992 } 00993 ns = 2*times + 1; 00994 nkv.Last() = new KnotVector(2, ns); 00995 KnotVector &lkv = *nkv.Last(); 00996 lkv[0] = lkv[1] = lkv[2] = 0.0; 00997 for (int i = 1; i < times; i++) 00998 lkv[2*i+1] = lkv[2*i+2] = i; 00999 lkv[ns] = lkv[ns+1] = lkv[ns+2] = times; 01000 lkv.GetElements(); 01001 NURBSPatch *newpatch = new NURBSPatch(nkv, 4); 01002 delete nkv.Last(); 01003 01004 DenseMatrix T(3), T2(3); 01005 Vector u(NULL, 3), v(NULL, 3); 01006 01007 NURBSPatch::Get3DRotationMatrix(n, ang, 1., T); 01008 double c = cos(ang/2); 01009 NURBSPatch::Get3DRotationMatrix(n, ang/2, 1./c, T2); 01010 T2 *= c; 01011 01012 double *op = patch.data, *np; 01013 for (int i = 0; i < size; i++) 01014 { 01015 np = newpatch->data + 4*i; 01016 for (int j = 0; j < 4; j++) 01017 np[j] = op[j]; 01018 for (int j = 0; j < times; j++) 01019 { 01020 u.SetData(np); 01021 v.SetData(np += 4*size); 01022 T2.Mult(u, v); 01023 v[3] = c*u[3]; 01024 v.SetData(np += 4*size); 01025 T.Mult(u, v); 01026 v[3] = u[3]; 01027 } 01028 op += 4; 01029 } 01030 01031 return newpatch; 01032 } 01033 01034 01035 // from mesh.cpp 01036 extern void skip_comment_lines(istream &, const char); 01037 01038 NURBSExtension::NURBSExtension(istream &input) 01039 { 01040 // Read topology 01041 patchTopo = new Mesh; 01042 patchTopo->LoadPatchTopo(input, edge_to_knot); 01043 own_topo = 1; 01044 01045 CheckPatches(); 01046 // CheckBdrPatches(); 01047 01048 skip_comment_lines(input, '#'); 01049 01050 // Read knotvectors or patches 01051 string ident; 01052 input >> ws >> ident; // 'knotvectors' or 'patches' 01053 if (ident == "knotvectors") 01054 { 01055 input >> NumOfKnotVectors; 01056 knotVectors.SetSize(NumOfKnotVectors); 01057 for (int i = 0; i < NumOfKnotVectors; i++) 01058 { 01059 knotVectors[i] = new KnotVector(input); 01060 if (knotVectors[i]->GetOrder() != knotVectors[0]->GetOrder()) 01061 mfem_error("NURBSExtension::NURBSExtension :\n" 01062 " Variable orders are not supported!"); 01063 } 01064 Order = knotVectors[0]->GetOrder(); 01065 } 01066 else if (ident == "patches") 01067 { 01068 patches.SetSize(GetNP()); 01069 for (int p = 0; p < patches.Size(); p++) 01070 { 01071 skip_comment_lines(input, '#'); 01072 patches[p] = new NURBSPatch(input); 01073 } 01074 01075 NumOfKnotVectors = 0; 01076 for (int i = 0; i < patchTopo->GetNEdges(); i++) 01077 if (NumOfKnotVectors < KnotInd(i)) 01078 NumOfKnotVectors = KnotInd(i); 01079 NumOfKnotVectors++; 01080 knotVectors.SetSize(NumOfKnotVectors); 01081 knotVectors = NULL; 01082 01083 Array<int> edges, oedge; 01084 for (int p = 0; p < patches.Size(); p++) 01085 { 01086 patchTopo->GetElementEdges(p, edges, oedge); 01087 if (Dimension() == 2) 01088 { 01089 if (knotVectors[KnotInd(edges[0])] == NULL) 01090 knotVectors[KnotInd(edges[0])] = 01091 new KnotVector(*patches[p]->GetKV(0)); 01092 if (knotVectors[KnotInd(edges[1])] == NULL) 01093 knotVectors[KnotInd(edges[1])] = 01094 new KnotVector(*patches[p]->GetKV(1)); 01095 } 01096 else 01097 { 01098 if (knotVectors[KnotInd(edges[0])] == NULL) 01099 knotVectors[KnotInd(edges[0])] = 01100 new KnotVector(*patches[p]->GetKV(0)); 01101 if (knotVectors[KnotInd(edges[3])] == NULL) 01102 knotVectors[KnotInd(edges[3])] = 01103 new KnotVector(*patches[p]->GetKV(1)); 01104 if (knotVectors[KnotInd(edges[8])] == NULL) 01105 knotVectors[KnotInd(edges[8])] = 01106 new KnotVector(*patches[p]->GetKV(2)); 01107 } 01108 } 01109 Order = knotVectors[0]->GetOrder(); 01110 } 01111 01112 GenerateOffsets(); 01113 CountElements(); 01114 CountBdrElements(); 01115 // NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs 01116 01117 skip_comment_lines(input, '#'); 01118 01119 // Check for a list of mesh elements 01120 if (patches.Size() == 0) 01121 input >> ws >> ident; 01122 if (patches.Size() == 0 && ident == "mesh_elements") 01123 { 01124 input >> NumOfActiveElems; 01125 activeElem.SetSize(GetGNE()); 01126 activeElem = false; 01127 int glob_elem; 01128 for (int i = 0; i < NumOfActiveElems; i++) 01129 { 01130 input >> glob_elem; 01131 activeElem[glob_elem] = true; 01132 } 01133 01134 skip_comment_lines(input, '#'); 01135 input >> ws >> ident; 01136 } 01137 else 01138 { 01139 NumOfActiveElems = NumOfElements; 01140 activeElem.SetSize(NumOfElements); 01141 activeElem = true; 01142 } 01143 01144 GenerateActiveVertices(); 01145 GenerateElementDofTable(); 01146 GenerateActiveBdrElems(); 01147 GenerateBdrElementDofTable(); 01148 01149 if (patches.Size() == 0) 01150 { 01151 // weights 01152 if (ident == "weights") 01153 { 01154 weights.Load(input, GetNDof()); 01155 } 01156 else // e.g. ident = "unitweights" or "autoweights" 01157 { 01158 weights.SetSize(GetNDof()); 01159 weights = 1.0; 01160 } 01161 } 01162 } 01163 01164 NURBSExtension::NURBSExtension(NURBSExtension *parent, int Order_) 01165 { 01166 Order = Order_; 01167 01168 patchTopo = parent->patchTopo; 01169 own_topo = 0; 01170 01171 parent->edge_to_knot.Copy(edge_to_knot); 01172 01173 NumOfKnotVectors = parent->GetNKV(); 01174 knotVectors.SetSize(NumOfKnotVectors); 01175 for (int i = 0; i < NumOfKnotVectors; i++) 01176 { 01177 knotVectors[i] = 01178 parent->GetKnotVector(i)->DegreeElevate(Order-parent->GetOrder()); 01179 } 01180 01181 // copy some data from parent 01182 NumOfElements = parent->NumOfElements; 01183 NumOfBdrElements = parent->NumOfBdrElements; 01184 01185 GenerateOffsets(); // dof offsets will be different from parent 01186 01187 NumOfActiveVertices = parent->NumOfActiveVertices; 01188 NumOfActiveElems = parent->NumOfActiveElems; 01189 NumOfActiveBdrElems = parent->NumOfActiveBdrElems; 01190 parent->activeVert.Copy(activeVert); 01191 parent->activeElem.Copy(activeElem); 01192 parent->activeBdrElem.Copy(activeBdrElem); 01193 01194 GenerateElementDofTable(); 01195 GenerateBdrElementDofTable(); 01196 01197 weights.SetSize(GetNDof()); 01198 weights = 1.0; 01199 } 01200 01201 NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces) 01202 { 01203 NURBSExtension *parent = mesh_array[0]->NURBSext; 01204 01205 if (!parent->own_topo) 01206 mfem_error("NURBSExtension::NURBSExtension :\n" 01207 " parent does not own the patch topology!"); 01208 patchTopo = parent->patchTopo; 01209 own_topo = 1; 01210 parent->own_topo = 0; 01211 01212 parent->edge_to_knot.Copy(edge_to_knot); 01213 01214 Order = parent->GetOrder(); 01215 01216 NumOfKnotVectors = parent->GetNKV(); 01217 knotVectors.SetSize(NumOfKnotVectors); 01218 for (int i = 0; i < NumOfKnotVectors; i++) 01219 { 01220 knotVectors[i] = new KnotVector(*parent->GetKnotVector(i)); 01221 } 01222 01223 GenerateOffsets(); 01224 CountElements(); 01225 CountBdrElements(); 01226 01227 // assuming the meshes define a partitioning of all the elements 01228 NumOfActiveElems = NumOfElements; 01229 activeElem.SetSize(NumOfElements); 01230 activeElem = true; 01231 01232 GenerateActiveVertices(); 01233 GenerateElementDofTable(); 01234 GenerateActiveBdrElems(); 01235 GenerateBdrElementDofTable(); 01236 01237 weights.SetSize(GetNDof()); 01238 MergeWeights(mesh_array, num_pieces); 01239 } 01240 01241 NURBSExtension::~NURBSExtension() 01242 { 01243 delete bel_dof; 01244 delete el_dof; 01245 01246 for (int i = 0; i < knotVectors.Size(); i++) 01247 delete knotVectors[i]; 01248 01249 for (int i = 0; i < patches.Size(); i++) 01250 delete patches[i]; 01251 01252 if (own_topo) 01253 delete patchTopo; 01254 } 01255 01256 void NURBSExtension::Print(ostream &out) const 01257 { 01258 patchTopo->PrintTopo(out, edge_to_knot); 01259 out << "\nknotvectors\n" << NumOfKnotVectors << '\n'; 01260 for (int i = 0; i < NumOfKnotVectors; i++) 01261 { 01262 knotVectors[i]->Print(out); 01263 } 01264 01265 if (NumOfActiveElems < NumOfElements) 01266 { 01267 out << "\nmesh_elements\n" << NumOfActiveElems << '\n'; 01268 for (int i = 0; i < NumOfElements; i++) 01269 if (activeElem[i]) 01270 out << i << '\n'; 01271 } 01272 01273 out << "\nweights\n"; 01274 weights.Print(out, 1); 01275 } 01276 01277 void NURBSExtension::PrintCharacteristics(ostream &out) 01278 { 01279 out << 01280 "NURBS Mesh entity sizes:\n" 01281 "Dimension = " << Dimension() << "\n" 01282 "Order = " << GetOrder() << "\n" 01283 "NumOfKnotVectors = " << GetNKV() << "\n" 01284 "NumOfPatches = " << GetNP() << "\n" 01285 "NumOfBdrPatches = " << GetNBP() << "\n" 01286 "NumOfVertices = " << GetGNV() << "\n" 01287 "NumOfElements = " << GetGNE() << "\n" 01288 "NumOfBdrElements = " << GetGNBE() << "\n" 01289 "NumOfDofs = " << GetNTotalDof() << "\n" 01290 "NumOfActiveVertices = " << GetNV() << "\n" 01291 "NumOfActiveElems = " << GetNE() << "\n" 01292 "NumOfActiveBdrElems = " << GetNBE() << "\n" 01293 "NumOfActiveDofs = " << GetNDof() << '\n'; 01294 for (int i = 0; i < NumOfKnotVectors; i++) 01295 { 01296 out << ' ' << i + 1 << ") "; 01297 knotVectors[i]->Print(out); 01298 } 01299 cout << endl; 01300 } 01301 01302 void NURBSExtension::GenerateActiveVertices() 01303 { 01304 int vert[8], nv, g_el, nx, ny, nz, dim = Dimension(); 01305 01306 NURBSPatchMap p2g(this); 01307 KnotVector *kv[3]; 01308 01309 g_el = 0; 01310 activeVert.SetSize(GetGNV()); 01311 activeVert = -1; 01312 for (int p = 0; p < GetNP(); p++) 01313 { 01314 p2g.SetPatchVertexMap(p, kv); 01315 01316 nx = p2g.nx(); 01317 ny = p2g.ny(); 01318 nz = (dim == 3) ? p2g.nz() : 1; 01319 01320 for (int k = 0; k < nz; k++) 01321 { 01322 for (int j = 0; j < ny; j++) 01323 { 01324 for (int i = 0; i < nx; i++) 01325 { 01326 if (activeElem[g_el]) 01327 { 01328 if (dim == 2) 01329 { 01330 vert[0] = p2g(i, j ); 01331 vert[1] = p2g(i+1,j ); 01332 vert[2] = p2g(i+1,j+1); 01333 vert[3] = p2g(i, j+1); 01334 nv = 4; 01335 } 01336 else 01337 { 01338 vert[0] = p2g(i, j, k); 01339 vert[1] = p2g(i+1,j, k); 01340 vert[2] = p2g(i+1,j+1,k); 01341 vert[3] = p2g(i, j+1,k); 01342 01343 vert[4] = p2g(i, j, k+1); 01344 vert[5] = p2g(i+1,j, k+1); 01345 vert[6] = p2g(i+1,j+1,k+1); 01346 vert[7] = p2g(i, j+1,k+1); 01347 nv = 8; 01348 } 01349 01350 for (int v = 0; v < nv; v++) 01351 activeVert[vert[v]] = 1; 01352 } 01353 g_el++; 01354 } 01355 } 01356 } 01357 } 01358 01359 NumOfActiveVertices = 0; 01360 for (int i = 0; i < GetGNV(); i++) 01361 if (activeVert[i] == 1) 01362 activeVert[i] = NumOfActiveVertices++; 01363 } 01364 01365 void NURBSExtension::GenerateActiveBdrElems() 01366 { 01367 int dim = Dimension(); 01368 Array<KnotVector *> kv(dim); 01369 01370 activeBdrElem.SetSize(GetGNBE()); 01371 if (GetGNE() == GetNE()) 01372 { 01373 activeBdrElem = true; 01374 NumOfActiveBdrElems = GetGNBE(); 01375 return; 01376 } 01377 activeBdrElem = false; 01378 NumOfActiveBdrElems = 0; 01379 // the mesh will generate the actual boundary including boundary 01380 // elements that are not on boundary patches. we use this for 01381 // visialization of processor boundaries 01382 01383 // TODO: generate actual boundary? 01384 } 01385 01386 01387 void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces) 01388 { 01389 Array<int> lelem_elem; 01390 01391 for (int i = 0; i < num_pieces; i++) 01392 { 01393 NURBSExtension *lext = mesh_array[i]->NURBSext; 01394 01395 lext->GetElementLocalToGlobal(lelem_elem); 01396 01397 for (int lel = 0; lel < lext->GetNE(); lel++) 01398 { 01399 int gel = lelem_elem[lel]; 01400 01401 int nd = el_dof->RowSize(gel); 01402 int *gdofs = el_dof->GetRow(gel); 01403 int *ldofs = lext->el_dof->GetRow(lel); 01404 for (int j = 0; j < nd; j++) 01405 weights(gdofs[j]) = lext->weights(ldofs[j]); 01406 } 01407 } 01408 } 01409 01410 void NURBSExtension::MergeGridFunctions( 01411 GridFunction *gf_array[], int num_pieces, GridFunction &merged) 01412 { 01413 FiniteElementSpace *gfes = merged.FESpace(); 01414 Array<int> lelem_elem, dofs; 01415 Vector lvec; 01416 01417 for (int i = 0; i < num_pieces; i++) 01418 { 01419 FiniteElementSpace *lfes = gf_array[i]->FESpace(); 01420 NURBSExtension *lext = lfes->GetMesh()->NURBSext; 01421 01422 lext->GetElementLocalToGlobal(lelem_elem); 01423 01424 for (int lel = 0; lel < lext->GetNE(); lel++) 01425 { 01426 lfes->GetElementVDofs(lel, dofs); 01427 gf_array[i]->GetSubVector(dofs, lvec); 01428 01429 gfes->GetElementVDofs(lelem_elem[lel], dofs); 01430 merged.SetSubVector(dofs, lvec); 01431 } 01432 } 01433 } 01434 01435 void NURBSExtension::CheckPatches() 01436 { 01437 Array<int> edges; 01438 Array<int> oedge; 01439 01440 for (int p = 0; p < GetNP(); p++) 01441 { 01442 patchTopo->GetElementEdges(p, edges, oedge); 01443 01444 for (int i = 0; i < edges.Size(); i++) 01445 { 01446 edges[i] = edge_to_knot[edges[i]]; 01447 if (oedge[i] < 0) 01448 edges[i] = -1 - edges[i]; 01449 } 01450 01451 if ((Dimension() == 2 && 01452 (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) || 01453 01454 (Dimension() == 3 && 01455 (edges[0] != edges[2] || edges[0] != edges[4] || 01456 edges[0] != edges[6] || edges[1] != edges[3] || 01457 edges[1] != edges[5] || edges[1] != edges[7] || 01458 edges[8] != edges[9] || edges[8] != edges[10] || 01459 edges[8] != edges[11]))) 01460 { 01461 cerr << "NURBSExtension::CheckPatch (patch = " << p 01462 << ")\n Inconsistent edge-to-knot mapping!\n"; 01463 mfem_error(); 01464 } 01465 01466 if ((Dimension() == 2 && 01467 (edges[0] < 0 || edges[1] < 0)) || 01468 01469 (Dimension() == 3 && 01470 (edges[0] < 0 || edges[3] < 0 || edges[8] < 0))) 01471 { 01472 cerr << "NURBSExtension::CheckPatch (patch = " << p 01473 << ") : Bad orientation!\n"; 01474 mfem_error(); 01475 } 01476 } 01477 } 01478 01479 void NURBSExtension::CheckBdrPatches() 01480 { 01481 Array<int> edges; 01482 Array<int> oedge; 01483 01484 for (int p = 0; p < GetNBP(); p++) 01485 { 01486 patchTopo->GetBdrElementEdges(p, edges, oedge); 01487 01488 for (int i = 0; i < edges.Size(); i++) 01489 { 01490 edges[i] = edge_to_knot[edges[i]]; 01491 if (oedge[i] < 0) 01492 edges[i] = -1 - edges[i]; 01493 } 01494 01495 if ((Dimension() == 2 && (edges[0] < 0)) || 01496 (Dimension() == 3 && (edges[0] < 0 || edges[1] < 0))) 01497 { 01498 cerr << "NURBSExtension::CheckBdrPatch (boundary patch = " 01499 << p << ") : Bad orientation!\n"; 01500 mfem_error(); 01501 } 01502 } 01503 } 01504 01505 void NURBSExtension::GetPatchKnotVectors(int p, Array<KnotVector *> &kv) 01506 { 01507 Array<int> edges; 01508 Array<int> orient; 01509 01510 kv.SetSize(Dimension()); 01511 patchTopo->GetElementEdges(p, edges, orient); 01512 01513 if (Dimension() == 2) 01514 { 01515 kv[0] = KnotVec(edges[0]); 01516 kv[1] = KnotVec(edges[1]); 01517 } 01518 else 01519 { 01520 kv[0] = KnotVec(edges[0]); 01521 kv[1] = KnotVec(edges[3]); 01522 kv[2] = KnotVec(edges[8]); 01523 } 01524 } 01525 01526 void NURBSExtension::GetBdrPatchKnotVectors(int p, Array<KnotVector *> &kv) 01527 { 01528 Array<int> edges; 01529 Array<int> orient; 01530 01531 kv.SetSize(Dimension() - 1); 01532 patchTopo->GetBdrElementEdges(p, edges, orient); 01533 01534 if (Dimension() == 2) 01535 { 01536 kv[0] = KnotVec(edges[0]); 01537 } 01538 else 01539 { 01540 kv[0] = KnotVec(edges[0]); 01541 kv[1] = KnotVec(edges[1]); 01542 } 01543 } 01544 01545 void NURBSExtension::GenerateOffsets() 01546 { 01547 int nv = patchTopo->GetNV(); 01548 int ne = patchTopo->GetNEdges(); 01549 int nf = patchTopo->GetNFaces(); 01550 int np = patchTopo->GetNE(); 01551 int meshCounter, spaceCounter, dim = Dimension(); 01552 01553 Array<int> edges; 01554 Array<int> orient; 01555 01556 v_meshOffsets.SetSize(nv); 01557 e_meshOffsets.SetSize(ne); 01558 f_meshOffsets.SetSize(nf); 01559 p_meshOffsets.SetSize(np); 01560 01561 v_spaceOffsets.SetSize(nv); 01562 e_spaceOffsets.SetSize(ne); 01563 f_spaceOffsets.SetSize(nf); 01564 p_spaceOffsets.SetSize(np); 01565 01566 // Get vertex offsets 01567 for (meshCounter = 0; meshCounter < nv; meshCounter++) 01568 { 01569 v_meshOffsets[meshCounter] = meshCounter; 01570 v_spaceOffsets[meshCounter] = meshCounter; 01571 } 01572 spaceCounter = meshCounter; 01573 01574 // Get edge offsets 01575 for (int e = 0; e < ne; e++) 01576 { 01577 e_meshOffsets[e] = meshCounter; 01578 e_spaceOffsets[e] = spaceCounter; 01579 meshCounter += KnotVec(e)->GetNE() - 1; 01580 spaceCounter += KnotVec(e)->GetNCP() - 2; 01581 } 01582 01583 // Get face offsets 01584 for (int f = 0; f < nf; f++) 01585 { 01586 f_meshOffsets[f] = meshCounter; 01587 f_spaceOffsets[f] = spaceCounter; 01588 01589 patchTopo->GetFaceEdges(f, edges, orient); 01590 01591 meshCounter += 01592 (KnotVec(edges[0])->GetNE() - 1) * 01593 (KnotVec(edges[1])->GetNE() - 1); 01594 spaceCounter += 01595 (KnotVec(edges[0])->GetNCP() - 2) * 01596 (KnotVec(edges[1])->GetNCP() - 2); 01597 } 01598 01599 // Get patch offsets 01600 for (int p = 0; p < np; p++) 01601 { 01602 p_meshOffsets[p] = meshCounter; 01603 p_spaceOffsets[p] = spaceCounter; 01604 01605 patchTopo->GetElementEdges(p, edges, orient); 01606 01607 if (dim == 2) 01608 { 01609 meshCounter += 01610 (KnotVec(edges[0])->GetNE() - 1) * 01611 (KnotVec(edges[1])->GetNE() - 1); 01612 spaceCounter += 01613 (KnotVec(edges[0])->GetNCP() - 2) * 01614 (KnotVec(edges[1])->GetNCP() - 2); 01615 } 01616 else 01617 { 01618 meshCounter += 01619 (KnotVec(edges[0])->GetNE() - 1) * 01620 (KnotVec(edges[3])->GetNE() - 1) * 01621 (KnotVec(edges[8])->GetNE() - 1); 01622 spaceCounter += 01623 (KnotVec(edges[0])->GetNCP() - 2) * 01624 (KnotVec(edges[3])->GetNCP() - 2) * 01625 (KnotVec(edges[8])->GetNCP() - 2); 01626 } 01627 } 01628 NumOfVertices = meshCounter; 01629 NumOfDofs = spaceCounter; 01630 } 01631 01632 void NURBSExtension::CountElements() 01633 { 01634 int dim = Dimension(); 01635 Array<KnotVector *> kv(dim); 01636 01637 NumOfElements = 0; 01638 for (int p = 0; p < GetNP(); p++) 01639 { 01640 GetPatchKnotVectors(p, kv); 01641 01642 int ne = kv[0]->GetNE(); 01643 for (int d = 1; d < dim; d++) 01644 ne *= kv[d]->GetNE(); 01645 01646 NumOfElements += ne; 01647 } 01648 } 01649 01650 void NURBSExtension::CountBdrElements() 01651 { 01652 int dim = Dimension() - 1; 01653 Array<KnotVector *> kv(dim); 01654 01655 NumOfBdrElements = 0; 01656 for (int p = 0; p < GetNBP(); p++) 01657 { 01658 GetBdrPatchKnotVectors(p, kv); 01659 01660 int ne = kv[0]->GetNE(); 01661 for (int d = 1; d < dim; d++) 01662 ne *= kv[d]->GetNE(); 01663 01664 NumOfBdrElements += ne; 01665 } 01666 } 01667 01668 void NURBSExtension::GetElementTopo(Array<Element *> &elements) 01669 { 01670 elements.SetSize(GetNE()); 01671 01672 if (Dimension() == 2) 01673 { 01674 Get2DElementTopo(elements); 01675 } 01676 else 01677 { 01678 Get3DElementTopo(elements); 01679 } 01680 } 01681 01682 void NURBSExtension::Get2DElementTopo(Array<Element *> &elements) 01683 { 01684 int el = 0; 01685 int eg = 0; 01686 int ind[4]; 01687 NURBSPatchMap p2g(this); 01688 KnotVector *kv[2]; 01689 01690 for (int p = 0; p < GetNP(); p++) 01691 { 01692 p2g.SetPatchVertexMap(p, kv); 01693 int nx = p2g.nx(); 01694 int ny = p2g.ny(); 01695 01696 int patch_attr = patchTopo->GetAttribute(p); 01697 01698 for (int j = 0; j < ny; j++) 01699 { 01700 for (int i = 0; i < nx; i++) 01701 { 01702 if (activeElem[eg]) 01703 { 01704 ind[0] = activeVert[p2g(i, j )]; 01705 ind[1] = activeVert[p2g(i+1,j )]; 01706 ind[2] = activeVert[p2g(i+1,j+1)]; 01707 ind[3] = activeVert[p2g(i, j+1)]; 01708 01709 elements[el] = new Quadrilateral(ind, patch_attr); 01710 el++; 01711 } 01712 eg++; 01713 } 01714 } 01715 } 01716 } 01717 01718 void NURBSExtension::Get3DElementTopo(Array<Element *> &elements) 01719 { 01720 int el = 0; 01721 int eg = 0; 01722 int ind[8]; 01723 NURBSPatchMap p2g(this); 01724 KnotVector *kv[3]; 01725 01726 for (int p = 0; p < GetNP(); p++) 01727 { 01728 p2g.SetPatchVertexMap(p, kv); 01729 int nx = p2g.nx(); 01730 int ny = p2g.ny(); 01731 int nz = p2g.nz(); 01732 01733 int patch_attr = patchTopo->GetAttribute(p); 01734 01735 for (int k = 0; k < nz; k++) 01736 { 01737 for (int j = 0; j < ny; j++) 01738 { 01739 for (int i = 0; i < nx; i++) 01740 { 01741 if (activeElem[eg]) 01742 { 01743 ind[0] = activeVert[p2g(i, j, k)]; 01744 ind[1] = activeVert[p2g(i+1,j, k)]; 01745 ind[2] = activeVert[p2g(i+1,j+1,k)]; 01746 ind[3] = activeVert[p2g(i, j+1,k)]; 01747 01748 ind[4] = activeVert[p2g(i, j, k+1)]; 01749 ind[5] = activeVert[p2g(i+1,j, k+1)]; 01750 ind[6] = activeVert[p2g(i+1,j+1,k+1)]; 01751 ind[7] = activeVert[p2g(i, j+1,k+1)]; 01752 01753 elements[el] = new Hexahedron(ind, patch_attr); 01754 el++; 01755 } 01756 eg++; 01757 } 01758 } 01759 } 01760 } 01761 } 01762 01763 void NURBSExtension::GetBdrElementTopo(Array<Element *> &boundary) 01764 { 01765 boundary.SetSize(GetNBE()); 01766 01767 if (Dimension() == 2) 01768 { 01769 Get2DBdrElementTopo(boundary); 01770 } 01771 else 01772 { 01773 Get3DBdrElementTopo(boundary); 01774 } 01775 } 01776 01777 void NURBSExtension::Get2DBdrElementTopo(Array<Element *> &boundary) 01778 { 01779 int g_be, l_be; 01780 int ind[2], okv[1]; 01781 NURBSPatchMap p2g(this); 01782 KnotVector *kv[1]; 01783 01784 g_be = l_be = 0; 01785 for (int b = 0; b < GetNBP(); b++) 01786 { 01787 p2g.SetBdrPatchVertexMap(b, kv, okv); 01788 int nx = p2g.nx(); 01789 01790 int bdr_patch_attr = patchTopo->GetBdrAttribute(b); 01791 01792 for (int i = 0; i < nx; i++) 01793 { 01794 if (activeBdrElem[g_be]) 01795 { 01796 int _i = (okv[0] >= 0) ? i : (nx - 1 - i); 01797 ind[0] = activeVert[p2g[_i ]]; 01798 ind[1] = activeVert[p2g[_i+1]]; 01799 01800 boundary[l_be] = new Segment(ind, bdr_patch_attr); 01801 l_be++; 01802 } 01803 g_be++; 01804 } 01805 } 01806 } 01807 01808 void NURBSExtension::Get3DBdrElementTopo(Array<Element *> &boundary) 01809 { 01810 int g_be, l_be; 01811 int ind[4], okv[2]; 01812 NURBSPatchMap p2g(this); 01813 KnotVector *kv[2]; 01814 01815 g_be = l_be = 0; 01816 for (int b = 0; b < GetNBP(); b++) 01817 { 01818 p2g.SetBdrPatchVertexMap(b, kv, okv); 01819 int nx = p2g.nx(); 01820 int ny = p2g.ny(); 01821 01822 int bdr_patch_attr = patchTopo->GetBdrAttribute(b); 01823 01824 for (int j = 0; j < ny; j++) 01825 { 01826 int _j = (okv[1] >= 0) ? j : (ny - 1 - j); 01827 for (int i = 0; i < nx; i++) 01828 { 01829 if (activeBdrElem[g_be]) 01830 { 01831 int _i = (okv[0] >= 0) ? i : (nx - 1 - i); 01832 ind[0] = activeVert[p2g(_i, _j )]; 01833 ind[1] = activeVert[p2g(_i+1,_j )]; 01834 ind[2] = activeVert[p2g(_i+1,_j+1)]; 01835 ind[3] = activeVert[p2g(_i, _j+1)]; 01836 01837 boundary[l_be] = new Quadrilateral(ind, bdr_patch_attr); 01838 l_be++; 01839 } 01840 g_be++; 01841 } 01842 } 01843 } 01844 } 01845 01846 void NURBSExtension::GenerateElementDofTable() 01847 { 01848 activeDof.SetSize(GetNTotalDof()); 01849 activeDof = 0; 01850 01851 if (Dimension() == 2) 01852 { 01853 Generate2DElementDofTable(); 01854 } 01855 else 01856 { 01857 Generate3DElementDofTable(); 01858 } 01859 01860 NumOfActiveDofs = 0; 01861 for (int d = 0; d < GetNTotalDof(); d++) 01862 if (activeDof[d]) 01863 { 01864 NumOfActiveDofs++; 01865 activeDof[d] = NumOfActiveDofs; 01866 } 01867 01868 int *dof = el_dof->GetJ(); 01869 int ndof = el_dof->Size_of_connections(); 01870 for (int i = 0; i < ndof; i++) 01871 { 01872 dof[i] = activeDof[dof[i]] - 1; 01873 } 01874 } 01875 01876 void NURBSExtension::Generate2DElementDofTable() 01877 { 01878 int el = 0; 01879 int eg = 0; 01880 KnotVector *kv[2]; 01881 NURBSPatchMap p2g(this); 01882 01883 el_dof = new Table(NumOfActiveElems, (Order + 1)*(Order + 1)); 01884 el_to_patch.SetSize(NumOfActiveElems); 01885 el_to_IJK.SetSize(NumOfActiveElems, 2); 01886 01887 for (int p = 0; p < GetNP(); p++) 01888 { 01889 p2g.SetPatchDofMap(p, kv); 01890 01891 // Load dofs 01892 for (int j = 0; j < kv[1]->GetNKS(); j++) 01893 { 01894 if (kv[1]->isElement(j)) 01895 { 01896 for (int i = 0; i < kv[0]->GetNKS(); i++) 01897 { 01898 if (kv[0]->isElement(i)) 01899 { 01900 if (activeElem[eg]) 01901 { 01902 int *dofs = el_dof->GetRow(el); 01903 int idx = 0; 01904 for (int jj = 0; jj <= Order; jj++) 01905 { 01906 for (int ii = 0; ii <= Order; ii++) 01907 { 01908 dofs[idx] = p2g(i+ii,j+jj); 01909 activeDof[dofs[idx]] = 1; 01910 idx++; 01911 } 01912 } 01913 el_to_patch[el] = p; 01914 el_to_IJK(el,0) = i; 01915 el_to_IJK(el,1) = j; 01916 01917 el++; 01918 } 01919 eg++; 01920 } 01921 } 01922 } 01923 } 01924 } 01925 } 01926 01927 void NURBSExtension::Generate3DElementDofTable() 01928 { 01929 int el = 0; 01930 int eg = 0; 01931 int idx; 01932 KnotVector *kv[3]; 01933 NURBSPatchMap p2g(this); 01934 01935 el_dof = new Table(NumOfActiveElems, (Order + 1)*(Order + 1)*(Order + 1)); 01936 el_to_patch.SetSize(NumOfActiveElems); 01937 el_to_IJK.SetSize(NumOfActiveElems, 3); 01938 01939 for (int p = 0; p < GetNP(); p++) 01940 { 01941 p2g.SetPatchDofMap(p, kv); 01942 01943 // Load dofs 01944 for (int k = 0; k < kv[2]->GetNKS(); k++) 01945 { 01946 if (kv[2]->isElement(k)) 01947 { 01948 for (int j = 0; j < kv[1]->GetNKS(); j++) 01949 { 01950 if (kv[1]->isElement(j)) 01951 { 01952 for (int i = 0; i < kv[0]->GetNKS(); i++) 01953 { 01954 if (kv[0]->isElement(i)) 01955 { 01956 if (activeElem[eg]) 01957 { 01958 idx = 0; 01959 int *dofs = el_dof->GetRow(el); 01960 for (int kk = 0; kk <= Order; kk++) 01961 { 01962 for (int jj = 0; jj <= Order; jj++) 01963 { 01964 for (int ii = 0; ii <= Order; ii++) 01965 { 01966 dofs[idx] = p2g(i+ii,j+jj,k+kk); 01967 activeDof[dofs[idx]] = 1; 01968 idx++; 01969 } 01970 } 01971 } 01972 01973 el_to_patch[el] = p; 01974 el_to_IJK(el,0) = i; 01975 el_to_IJK(el,1) = j; 01976 el_to_IJK(el,2) = k; 01977 01978 el++; 01979 } 01980 eg++; 01981 } 01982 } 01983 } 01984 } 01985 } 01986 } 01987 } 01988 } 01989 01990 void NURBSExtension::GenerateBdrElementDofTable() 01991 { 01992 if (Dimension() == 2) 01993 { 01994 Generate2DBdrElementDofTable(); 01995 } 01996 else 01997 { 01998 Generate3DBdrElementDofTable(); 01999 } 02000 02001 int *dof = bel_dof->GetJ(); 02002 int ndof = bel_dof->Size_of_connections(); 02003 for (int i = 0; i < ndof; i++) 02004 { 02005 dof[i] = activeDof[dof[i]] - 1; 02006 } 02007 } 02008 02009 void NURBSExtension::Generate2DBdrElementDofTable() 02010 { 02011 int gbe = 0; 02012 int lbe = 0, okv[1]; 02013 KnotVector *kv[1]; 02014 NURBSPatchMap p2g(this); 02015 02016 bel_dof = new Table(NumOfActiveBdrElems, Order + 1); 02017 bel_to_patch.SetSize(NumOfActiveBdrElems); 02018 bel_to_IJK.SetSize(NumOfActiveBdrElems, 1); 02019 02020 for (int b = 0; b < GetNBP(); b++) 02021 { 02022 p2g.SetBdrPatchDofMap(b, kv, okv); 02023 int nx = p2g.nx(); // NCP-1 02024 02025 // Load dofs 02026 int nks0 = kv[0]->GetNKS(); 02027 for (int i = 0; i < nks0; i++) 02028 { 02029 if (kv[0]->isElement(i)) 02030 { 02031 if (activeBdrElem[gbe]) 02032 { 02033 int *dofs = bel_dof->GetRow(lbe); 02034 int idx = 0; 02035 for (int ii = 0; ii <= Order; ii++) 02036 { 02037 dofs[idx] = p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)]; 02038 idx++; 02039 } 02040 bel_to_patch[lbe] = b; 02041 bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i); 02042 lbe++; 02043 } 02044 gbe++; 02045 } 02046 } 02047 } 02048 } 02049 02050 void NURBSExtension::Generate3DBdrElementDofTable() 02051 { 02052 int gbe = 0; 02053 int lbe = 0, okv[2]; 02054 KnotVector *kv[2]; 02055 NURBSPatchMap p2g(this); 02056 02057 bel_dof = new Table(NumOfActiveBdrElems, (Order + 1)*(Order + 1)); 02058 bel_to_patch.SetSize(NumOfActiveBdrElems); 02059 bel_to_IJK.SetSize(NumOfActiveBdrElems, 2); 02060 02061 for (int b = 0; b < GetNBP(); b++) 02062 { 02063 p2g.SetBdrPatchDofMap(b, kv, okv); 02064 int nx = p2g.nx(); // NCP0-1 02065 int ny = p2g.ny(); // NCP1-1 02066 02067 // Load dofs 02068 int nks0 = kv[0]->GetNKS(); 02069 int nks1 = kv[1]->GetNKS(); 02070 for (int j = 0; j < nks1; j++) 02071 { 02072 if (kv[1]->isElement(j)) 02073 { 02074 for (int i = 0; i < nks0; i++) 02075 { 02076 if (kv[0]->isElement(i)) 02077 { 02078 if (activeBdrElem[gbe]) 02079 { 02080 int *dofs = bel_dof->GetRow(lbe); 02081 int idx = 0; 02082 for (int jj = 0; jj <= Order; jj++) 02083 { 02084 int _jj = (okv[1] >= 0) ? (j+jj) : (ny-j-jj); 02085 for (int ii = 0; ii <= Order; ii++) 02086 { 02087 int _ii = (okv[0] >= 0) ? (i+ii) : (nx-i-ii); 02088 dofs[idx] = p2g(_ii,_jj); 02089 idx++; 02090 } 02091 } 02092 bel_to_patch[lbe] = b; 02093 bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i); 02094 bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j); 02095 lbe++; 02096 } 02097 gbe++; 02098 } 02099 } 02100 } 02101 } 02102 } 02103 } 02104 02105 void NURBSExtension::GetVertexLocalToGlobal(Array<int> &lvert_vert) 02106 { 02107 lvert_vert.SetSize(GetNV()); 02108 for (int gv = 0; gv < GetGNV(); gv++) 02109 if (activeVert[gv] >= 0) 02110 lvert_vert[activeVert[gv]] = gv; 02111 } 02112 02113 void NURBSExtension::GetElementLocalToGlobal(Array<int> &lelem_elem) 02114 { 02115 lelem_elem.SetSize(GetNE()); 02116 for (int le = 0, ge = 0; ge < GetGNE(); ge++) 02117 if (activeElem[ge]) 02118 lelem_elem[le++] = ge; 02119 } 02120 02121 void NURBSExtension::LoadFE(int i, const FiniteElement *FE) 02122 { 02123 const NURBSFiniteElement *NURBSFE = 02124 dynamic_cast<const NURBSFiniteElement *>(FE); 02125 02126 if (NURBSFE->GetElement() != i) 02127 { 02128 Array<int> dofs; 02129 NURBSFE->SetIJK(el_to_IJK.GetRow(i)); 02130 if (el_to_patch[i] != NURBSFE->GetPatch()) 02131 { 02132 GetPatchKnotVectors(el_to_patch[i], NURBSFE->KnotVectors()); 02133 NURBSFE->SetPatch(el_to_patch[i]); 02134 } 02135 el_dof->GetRow(i, dofs); 02136 weights.GetSubVector(dofs, NURBSFE->Weights()); 02137 NURBSFE->SetElement(i); 02138 } 02139 } 02140 02141 void NURBSExtension::LoadBE(int i, const FiniteElement *BE) 02142 { 02143 const NURBSFiniteElement *NURBSFE = 02144 dynamic_cast<const NURBSFiniteElement *>(BE); 02145 02146 if (NURBSFE->GetElement() != i) 02147 { 02148 Array<int> dofs; 02149 NURBSFE->SetIJK(bel_to_IJK.GetRow(i)); 02150 if (bel_to_patch[i] != NURBSFE->GetPatch()) 02151 { 02152 GetBdrPatchKnotVectors(bel_to_patch[i], NURBSFE->KnotVectors()); 02153 NURBSFE->SetPatch(bel_to_patch[i]); 02154 } 02155 bel_dof->GetRow(i, dofs); 02156 weights.GetSubVector(dofs, NURBSFE->Weights()); 02157 NURBSFE->SetElement(i); 02158 } 02159 } 02160 02161 void NURBSExtension::ConvertToPatches(const Vector &Nodes) 02162 { 02163 delete el_dof; 02164 delete bel_dof; 02165 02166 if (patches.Size() == 0) 02167 GetPatchNets(Nodes); 02168 } 02169 02170 void NURBSExtension::SetCoordsFromPatches(Vector &Nodes) 02171 { 02172 if (patches.Size() == 0) return; 02173 02174 SetSolutionVector(Nodes); 02175 patches.SetSize(0); 02176 } 02177 02178 void NURBSExtension::SetKnotsFromPatches() 02179 { 02180 if (patches.Size() == 0) 02181 mfem_error("NURBSExtension::SetKnotsFromPatches :" 02182 " No patches available!"); 02183 02184 Array<KnotVector *> kv; 02185 02186 for (int p = 0; p < patches.Size(); p++) 02187 { 02188 GetPatchKnotVectors(p, kv); 02189 02190 for (int i = 0; i < kv.Size(); i++) 02191 *kv[i] = *patches[p]->GetKV(i); 02192 } 02193 02194 Order = knotVectors[0]->GetOrder(); 02195 for (int i = 0; i < NumOfKnotVectors; i++) 02196 { 02197 if (Order != knotVectors[i]->GetOrder()) 02198 mfem_error("NURBSExtension::Reset :\n" 02199 " Variable orders are not supported!"); 02200 } 02201 02202 GenerateOffsets(); 02203 CountElements(); 02204 CountBdrElements(); 02205 02206 // all elements must be active 02207 NumOfActiveElems = NumOfElements; 02208 activeElem.SetSize(NumOfElements); 02209 activeElem = true; 02210 02211 GenerateActiveVertices(); 02212 GenerateElementDofTable(); 02213 GenerateActiveBdrElems(); 02214 GenerateBdrElementDofTable(); 02215 } 02216 02217 void NURBSExtension::DegreeElevate(int t) 02218 { 02219 for (int p = 0; p < patches.Size(); p++) 02220 { 02221 patches[p]->DegreeElevate(t); 02222 } 02223 } 02224 02225 void NURBSExtension::UniformRefinement() 02226 { 02227 for (int p = 0; p < patches.Size(); p++) 02228 { 02229 patches[p]->UniformRefinement(); 02230 } 02231 } 02232 02233 void NURBSExtension::KnotInsert(Array<KnotVector *> &kv) 02234 { 02235 Array<int> edges; 02236 Array<int> orient; 02237 02238 Array<KnotVector *> pkv(Dimension()); 02239 02240 for (int p = 0; p < patches.Size(); p++) 02241 { 02242 patchTopo->GetElementEdges(p, edges, orient); 02243 02244 if (Dimension()==2) 02245 { 02246 pkv[0] = kv[KnotInd(edges[0])]; 02247 pkv[1] = kv[KnotInd(edges[1])]; 02248 } 02249 else 02250 { 02251 pkv[0] = kv[KnotInd(edges[0])]; 02252 pkv[1] = kv[KnotInd(edges[3])]; 02253 pkv[2] = kv[KnotInd(edges[8])]; 02254 } 02255 02256 patches[p]->KnotInsert(pkv); 02257 } 02258 } 02259 02260 void NURBSExtension::GetPatchNets(const Vector &coords) 02261 { 02262 if (Dimension() == 2) 02263 { 02264 Get2DPatchNets(coords); 02265 } 02266 else 02267 { 02268 Get3DPatchNets(coords); 02269 } 02270 } 02271 02272 void NURBSExtension::Get2DPatchNets(const Vector &coords) 02273 { 02274 Array<KnotVector *> kv(2); 02275 NURBSPatchMap p2g(this); 02276 02277 patches.SetSize(GetNP()); 02278 for (int p = 0; p < GetNP(); p++) 02279 { 02280 p2g.SetPatchDofMap(p, kv); 02281 patches[p] = new NURBSPatch(kv, 2+1); 02282 NURBSPatch &Patch = *patches[p]; 02283 02284 for (int j = 0; j < kv[1]->GetNCP(); j++) 02285 { 02286 for (int i = 0; i < kv[0]->GetNCP(); i++) 02287 { 02288 int l = p2g(i,j); 02289 for (int d = 0; d < 2; d++) 02290 { 02291 Patch(i,j,d) = coords(l*2 + d)*weights(l); 02292 } 02293 Patch(i,j,2) = weights(l); 02294 } 02295 } 02296 } 02297 } 02298 02299 void NURBSExtension::Get3DPatchNets(const Vector &coords) 02300 { 02301 Array<KnotVector *> kv(3); 02302 NURBSPatchMap p2g(this); 02303 02304 patches.SetSize(GetNP()); 02305 for (int p = 0; p < GetNP(); p++) 02306 { 02307 p2g.SetPatchDofMap(p, kv); 02308 patches[p] = new NURBSPatch(kv, 3+1); 02309 NURBSPatch &Patch = *patches[p]; 02310 02311 for (int k = 0; k < kv[2]->GetNCP(); k++) 02312 { 02313 for (int j = 0; j < kv[1]->GetNCP(); j++) 02314 { 02315 for (int i = 0; i < kv[0]->GetNCP(); i++) 02316 { 02317 int l = p2g(i,j,k); 02318 for (int d = 0; d < 3; d++) 02319 { 02320 Patch(i,j,k,d) = coords(l*3 + d)*weights(l); 02321 } 02322 Patch(i,j,k,3) = weights(l); 02323 } 02324 } 02325 } 02326 } 02327 } 02328 02329 void NURBSExtension::SetSolutionVector(Vector &coords) 02330 { 02331 if (Dimension() == 2) 02332 { 02333 Set2DSolutionVector(coords); 02334 } 02335 else 02336 { 02337 Set3DSolutionVector(coords); 02338 } 02339 } 02340 02341 void NURBSExtension::Set2DSolutionVector(Vector &coords) 02342 { 02343 Array<KnotVector *> kv(2); 02344 NURBSPatchMap p2g(this); 02345 02346 weights.SetSize(GetNDof()); 02347 for (int p = 0; p < GetNP(); p++) 02348 { 02349 p2g.SetPatchDofMap(p, kv); 02350 NURBSPatch &Patch = *patches[p]; 02351 02352 for (int j = 0; j < kv[1]->GetNCP(); j++) 02353 { 02354 for (int i = 0; i < kv[0]->GetNCP(); i++) 02355 { 02356 int l = p2g(i,j); 02357 for (int d = 0; d < 2; d++) 02358 { 02359 coords(l*2 + d) = Patch(i,j,d)/Patch(i,j,2); 02360 } 02361 weights(l) = Patch(i,j,2); 02362 } 02363 } 02364 delete patches[p]; 02365 } 02366 } 02367 02368 void NURBSExtension::Set3DSolutionVector(Vector &coords) 02369 { 02370 Array<KnotVector *> kv(3); 02371 NURBSPatchMap p2g(this); 02372 02373 weights.SetSize(GetNDof()); 02374 for (int p = 0; p < GetNP(); p++) 02375 { 02376 p2g.SetPatchDofMap(p, kv); 02377 NURBSPatch &Patch = *patches[p]; 02378 02379 for (int k = 0; k < kv[2]->GetNCP(); k++) 02380 { 02381 for (int j = 0; j < kv[1]->GetNCP(); j++) 02382 { 02383 for (int i = 0; i < kv[0]->GetNCP(); i++) 02384 { 02385 int l = p2g(i,j,k); 02386 for (int d = 0; d < 3; d++) 02387 { 02388 coords(l*3 + d) = Patch(i,j,k,d)/Patch(i,j,k,3); 02389 } 02390 weights(l) = Patch(i,j,k,3); 02391 } 02392 } 02393 } 02394 delete patches[p]; 02395 } 02396 } 02397 02398 02399 #ifdef MFEM_USE_MPI 02400 ParNURBSExtension::ParNURBSExtension(MPI_Comm comm, NURBSExtension *parent, 02401 int *part, const Array<bool> &active_bel) 02402 : gtopo(comm) 02403 { 02404 if (parent->NumOfActiveElems < parent->NumOfElements) 02405 // SetActive (BuildGroups?) and the way the weights are copied 02406 // do not support this case 02407 mfem_error("ParNURBSExtension::ParNURBSExtension :\n" 02408 " all elements in the parent must be active!"); 02409 02410 patchTopo = parent->patchTopo; 02411 // steal ownership of patchTopo from the 'parent' NURBS extension 02412 if (!parent->own_topo) 02413 mfem_error("ParNURBSExtension::ParNURBSExtension :\n" 02414 " parent does not own the patch topology!"); 02415 own_topo = 1; 02416 parent->own_topo = 0; 02417 02418 parent->edge_to_knot.Copy(edge_to_knot); 02419 02420 Order = parent->GetOrder(); 02421 02422 NumOfKnotVectors = parent->GetNKV(); 02423 knotVectors.SetSize(NumOfKnotVectors); 02424 for (int i = 0; i < NumOfKnotVectors; i++) 02425 { 02426 knotVectors[i] = new KnotVector(*parent->GetKnotVector(i)); 02427 } 02428 02429 GenerateOffsets(); 02430 CountElements(); 02431 CountBdrElements(); 02432 02433 partitioning = part; 02434 SetActive(partitioning, active_bel); 02435 02436 GenerateActiveVertices(); 02437 GenerateElementDofTable(); 02438 // GenerateActiveBdrElems(); // done by SetActive for now 02439 GenerateBdrElementDofTable(); 02440 02441 Table *serial_elem_dof = parent->GetElementDofTable(); 02442 BuildGroups(partitioning, *serial_elem_dof); 02443 02444 weights.SetSize(GetNDof()); 02445 // copy weights from parent 02446 for (int gel = 0, lel = 0; gel < GetGNE(); gel++) 02447 { 02448 if (activeElem[gel]) 02449 { 02450 int ndofs = el_dof->RowSize(lel); 02451 int *ldofs = el_dof->GetRow(lel); 02452 int *gdofs = serial_elem_dof->GetRow(gel); 02453 for (int i = 0; i < ndofs; i++) 02454 weights(ldofs[i]) = parent->weights(gdofs[i]); 02455 lel++; 02456 } 02457 } 02458 } 02459 02460 ParNURBSExtension::ParNURBSExtension(NURBSExtension *parent, 02461 ParNURBSExtension *par_parent) 02462 : gtopo(par_parent->gtopo.GetComm()) 02463 { 02464 // steal all data from parent 02465 Order = parent->Order; 02466 02467 patchTopo = parent->patchTopo; 02468 own_topo = parent->own_topo; 02469 parent->own_topo = 0; 02470 02471 Swap(edge_to_knot, parent->edge_to_knot); 02472 02473 NumOfKnotVectors = parent->NumOfKnotVectors; 02474 Swap(knotVectors, parent->knotVectors); 02475 02476 NumOfVertices = parent->NumOfVertices; 02477 NumOfElements = parent->NumOfElements; 02478 NumOfBdrElements = parent->NumOfBdrElements; 02479 NumOfDofs = parent->NumOfDofs; 02480 02481 Swap(v_meshOffsets, parent->v_meshOffsets); 02482 Swap(e_meshOffsets, parent->e_meshOffsets); 02483 Swap(f_meshOffsets, parent->f_meshOffsets); 02484 Swap(p_meshOffsets, parent->p_meshOffsets); 02485 02486 Swap(v_spaceOffsets, parent->v_spaceOffsets); 02487 Swap(e_spaceOffsets, parent->e_spaceOffsets); 02488 Swap(f_spaceOffsets, parent->f_spaceOffsets); 02489 Swap(p_spaceOffsets, parent->p_spaceOffsets); 02490 02491 NumOfActiveVertices = parent->NumOfActiveVertices; 02492 NumOfActiveElems = parent->NumOfActiveElems; 02493 NumOfActiveBdrElems = parent->NumOfActiveBdrElems; 02494 NumOfActiveDofs = parent->NumOfActiveDofs; 02495 02496 Swap(activeVert, parent->activeVert); 02497 Swap(activeElem, parent->activeElem); 02498 Swap(activeBdrElem, parent->activeBdrElem); 02499 Swap(activeDof, parent->activeDof); 02500 02501 el_dof = parent->el_dof; 02502 bel_dof = parent->bel_dof; 02503 parent->el_dof = parent->bel_dof = NULL; 02504 02505 Swap(el_to_patch, parent->el_to_patch); 02506 Swap(bel_to_patch, parent->bel_to_patch); 02507 Swap(el_to_IJK, parent->el_to_IJK); 02508 Swap(bel_to_IJK, parent->bel_to_IJK); 02509 02510 swap(&weights, &parent->weights); 02511 02512 delete parent; 02513 02514 partitioning = NULL; 02515 02516 #ifdef MFEM_DEBUG 02517 if (par_parent->partitioning == NULL) 02518 mfem_error("ParNURBSExtension::ParNURBSExtension :\n" 02519 " parent ParNURBSExtension has no partitioning!"); 02520 #endif 02521 02522 Table *serial_elem_dof = GetGlobalElementDofTable(); 02523 BuildGroups(par_parent->partitioning, *serial_elem_dof); 02524 delete serial_elem_dof; 02525 } 02526 02527 Table *ParNURBSExtension::GetGlobalElementDofTable() 02528 { 02529 if (Dimension() == 2) 02530 return Get2DGlobalElementDofTable(); 02531 else 02532 return Get3DGlobalElementDofTable(); 02533 } 02534 02535 Table *ParNURBSExtension::Get2DGlobalElementDofTable() 02536 { 02537 int el = 0; 02538 KnotVector *kv[2]; 02539 NURBSPatchMap p2g(this); 02540 02541 Table *gel_dof = new Table(GetGNE(), (Order + 1)*(Order + 1)); 02542 02543 for (int p = 0; p < GetNP(); p++) 02544 { 02545 p2g.SetPatchDofMap(p, kv); 02546 02547 // Load dofs 02548 for (int j = 0; j < kv[1]->GetNKS(); j++) 02549 { 02550 if (kv[1]->isElement(j)) 02551 { 02552 for (int i = 0; i < kv[0]->GetNKS(); i++) 02553 { 02554 if (kv[0]->isElement(i)) 02555 { 02556 int *dofs = gel_dof->GetRow(el); 02557 int idx = 0; 02558 for (int jj = 0; jj <= Order; jj++) 02559 { 02560 for (int ii = 0; ii <= Order; ii++) 02561 { 02562 dofs[idx] = p2g(i+ii,j+jj); 02563 idx++; 02564 } 02565 } 02566 el++; 02567 } 02568 } 02569 } 02570 } 02571 } 02572 return gel_dof; 02573 } 02574 02575 Table *ParNURBSExtension::Get3DGlobalElementDofTable() 02576 { 02577 int el = 0; 02578 KnotVector *kv[3]; 02579 NURBSPatchMap p2g(this); 02580 02581 Table *gel_dof = new Table(GetGNE(), (Order + 1)*(Order + 1)*(Order + 1)); 02582 02583 for (int p = 0; p < GetNP(); p++) 02584 { 02585 p2g.SetPatchDofMap(p, kv); 02586 02587 // Load dofs 02588 for (int k = 0; k < kv[2]->GetNKS(); k++) 02589 { 02590 if (kv[2]->isElement(k)) 02591 { 02592 for (int j = 0; j < kv[1]->GetNKS(); j++) 02593 { 02594 if (kv[1]->isElement(j)) 02595 { 02596 for (int i = 0; i < kv[0]->GetNKS(); i++) 02597 { 02598 if (kv[0]->isElement(i)) 02599 { 02600 int *dofs = gel_dof->GetRow(el); 02601 int idx = 0; 02602 for (int kk = 0; kk <= Order; kk++) 02603 { 02604 for (int jj = 0; jj <= Order; jj++) 02605 { 02606 for (int ii = 0; ii <= Order; ii++) 02607 { 02608 dofs[idx] = p2g(i+ii,j+jj,k+kk); 02609 idx++; 02610 } 02611 } 02612 } 02613 el++; 02614 } 02615 } 02616 } 02617 } 02618 } 02619 } 02620 } 02621 return gel_dof; 02622 } 02623 02624 void ParNURBSExtension::SetActive(int *_partitioning, 02625 const Array<bool> &active_bel) 02626 { 02627 activeElem.SetSize(GetGNE()); 02628 activeElem = false; 02629 NumOfActiveElems = 0; 02630 int MyRank = gtopo.MyRank(); 02631 for (int i = 0; i < GetGNE(); i++) 02632 if (_partitioning[i] == MyRank) 02633 { 02634 activeElem[i] = true; 02635 NumOfActiveElems++; 02636 } 02637 02638 active_bel.Copy(activeBdrElem); 02639 NumOfActiveBdrElems = 0; 02640 for (int i = 0; i < GetGNBE(); i++) 02641 if (activeBdrElem[i]) 02642 NumOfActiveBdrElems++; 02643 } 02644 02645 void ParNURBSExtension::BuildGroups(int *_partitioning, const Table &elem_dof) 02646 { 02647 Table dof_proc; 02648 02649 ListOfIntegerSets groups; 02650 IntegerSet group; 02651 02652 Transpose(elem_dof, dof_proc); // dof_proc is dof_elem 02653 // convert elements to processors 02654 for (int i = 0; i < dof_proc.Size_of_connections(); i++) 02655 dof_proc.GetJ()[i] = _partitioning[dof_proc.GetJ()[i]]; 02656 02657 // the first group is the local one 02658 int MyRank = gtopo.MyRank(); 02659 group.Recreate(1, &MyRank); 02660 groups.Insert(group); 02661 02662 int dof = 0; 02663 ldof_group.SetSize(GetNDof()); 02664 for (int d = 0; d < GetNTotalDof(); d++) 02665 if (activeDof[d]) 02666 { 02667 group.Recreate(dof_proc.RowSize(d), dof_proc.GetRow(d)); 02668 ldof_group[dof] = groups.Insert(group); 02669 02670 dof++; 02671 } 02672 02673 gtopo.Create(groups, 1822); 02674 } 02675 #endif 02676 02677 02678 void NURBSPatchMap::GetPatchKnotVectors(int p, KnotVector *kv[]) 02679 { 02680 Ext->patchTopo->GetElementVertices(p, verts); 02681 Ext->patchTopo->GetElementEdges(p, edges, oedge); 02682 if (Ext->Dimension() == 2) 02683 { 02684 kv[0] = Ext->KnotVec(edges[0]); 02685 kv[1] = Ext->KnotVec(edges[1]); 02686 } 02687 else 02688 { 02689 Ext->patchTopo->GetElementFaces(p, faces, oface); 02690 02691 kv[0] = Ext->KnotVec(edges[0]); 02692 kv[1] = Ext->KnotVec(edges[3]); 02693 kv[2] = Ext->KnotVec(edges[8]); 02694 } 02695 opatch = 0; 02696 } 02697 02698 void NURBSPatchMap::GetBdrPatchKnotVectors(int p, KnotVector *kv[], int *okv) 02699 { 02700 Ext->patchTopo->GetBdrElementVertices(p, verts); 02701 Ext->patchTopo->GetBdrElementEdges(p, edges, oedge); 02702 kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]); 02703 if (Ext->Dimension() == 3) 02704 { 02705 faces.SetSize(1); 02706 Ext->patchTopo->GetBdrElementFace(p, &faces[0], &opatch); 02707 kv[1] = Ext->KnotVec(edges[1], oedge[1], &okv[1]); 02708 } 02709 else 02710 opatch = oedge[0]; 02711 } 02712 02713 void NURBSPatchMap::SetPatchVertexMap(int p, KnotVector *kv[]) 02714 { 02715 GetPatchKnotVectors(p, kv); 02716 02717 I = kv[0]->GetNE() - 1; 02718 J = kv[1]->GetNE() - 1; 02719 02720 for (int i = 0; i < verts.Size(); i++) 02721 verts[i] = Ext->v_meshOffsets[verts[i]]; 02722 02723 for (int i = 0; i < edges.Size(); i++) 02724 edges[i] = Ext->e_meshOffsets[edges[i]]; 02725 02726 if (Ext->Dimension() == 3) 02727 { 02728 K = kv[2]->GetNE() - 1; 02729 02730 for (int i = 0; i < faces.Size(); i++) 02731 faces[i] = Ext->f_meshOffsets[faces[i]]; 02732 } 02733 02734 pOffset = Ext->p_meshOffsets[p]; 02735 } 02736 02737 void NURBSPatchMap::SetPatchDofMap(int p, KnotVector *kv[]) 02738 { 02739 GetPatchKnotVectors(p, kv); 02740 02741 I = kv[0]->GetNCP() - 2; 02742 J = kv[1]->GetNCP() - 2; 02743 02744 for (int i = 0; i < verts.Size(); i++) 02745 verts[i] = Ext->v_spaceOffsets[verts[i]]; 02746 02747 for (int i = 0; i < edges.Size(); i++) 02748 edges[i] = Ext->e_spaceOffsets[edges[i]]; 02749 02750 if (Ext->Dimension() == 3) 02751 { 02752 K = kv[2]->GetNCP() - 2; 02753 02754 for (int i = 0; i < faces.Size(); i++) 02755 faces[i] = Ext->f_spaceOffsets[faces[i]]; 02756 } 02757 02758 pOffset = Ext->p_spaceOffsets[p]; 02759 } 02760 02761 void NURBSPatchMap::SetBdrPatchVertexMap(int p, KnotVector *kv[], int *okv) 02762 { 02763 GetBdrPatchKnotVectors(p, kv, okv); 02764 02765 I = kv[0]->GetNE() - 1; 02766 02767 for (int i = 0; i < verts.Size(); i++) 02768 verts[i] = Ext->v_meshOffsets[verts[i]]; 02769 02770 if (Ext->Dimension() == 2) 02771 { 02772 pOffset = Ext->e_meshOffsets[edges[0]]; 02773 } 02774 else 02775 { 02776 J = kv[1]->GetNE() - 1; 02777 02778 for (int i = 0; i < edges.Size(); i++) 02779 edges[i] = Ext->e_meshOffsets[edges[i]]; 02780 02781 pOffset = Ext->f_meshOffsets[faces[0]]; 02782 } 02783 } 02784 02785 void NURBSPatchMap::SetBdrPatchDofMap(int p, KnotVector *kv[], int *okv) 02786 { 02787 GetBdrPatchKnotVectors(p, kv, okv); 02788 02789 I = kv[0]->GetNCP() - 2; 02790 02791 for (int i = 0; i < verts.Size(); i++) 02792 verts[i] = Ext->v_spaceOffsets[verts[i]]; 02793 02794 if (Ext->Dimension() == 2) 02795 { 02796 pOffset = Ext->e_spaceOffsets[edges[0]]; 02797 } 02798 else 02799 { 02800 J = kv[1]->GetNCP() - 2; 02801 02802 for (int i = 0; i < edges.Size(); i++) 02803 edges[i] = Ext->e_spaceOffsets[edges[i]]; 02804 02805 pOffset = Ext->f_spaceOffsets[faces[0]]; 02806 } 02807 }