MFEM v2.0
nurbs.cpp
Go to the documentation of this file.
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 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines