00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceBernstein.hpp"
00032 #include "SundanceADReal.hpp"
00033 #include "PlayaExceptions.hpp"
00034 #include "SundanceSpatialDerivSpecifier.hpp"
00035 #include "SundancePoint.hpp"
00036 #include "SundanceObjectWithVerbosity.hpp"
00037 #include "SundanceOut.hpp"
00038
00039 using namespace Sundance;
00040 using namespace Teuchos;
00041
00042
00043 Bernstein::Bernstein(int order)
00044 : order_(order)
00045 {
00046 TEUCHOS_TEST_FOR_EXCEPTION(order < 0, std::runtime_error,
00047 "invalid polynomial order=" << order
00048 << " in Bernstein ctor");
00049 }
00050
00051 bool Bernstein::supportsCellTypePair(
00052 const CellType& maximalCellType,
00053 const CellType& cellType
00054 ) const
00055 {
00056 switch(maximalCellType)
00057 {
00058 case LineCell:
00059 switch(cellType)
00060 {
00061 case LineCell:
00062 case PointCell:
00063 return true;
00064 default:
00065 return false;
00066 }
00067 case TriangleCell:
00068 switch(cellType)
00069 {
00070 case TriangleCell:
00071 case LineCell:
00072 case PointCell:
00073 return true;
00074 default:
00075 return false;
00076 }
00077 case TetCell:
00078 switch(cellType)
00079 {
00080 case TetCell:
00081 case TriangleCell:
00082 case LineCell:
00083 case PointCell:
00084 return true;
00085 default:
00086 return false;
00087 }
00088 default:
00089 return false;
00090 }
00091 }
00092
00093 void Bernstein::print(std::ostream& os) const
00094 {
00095 os << "Bernstein(" << order_ << ")";
00096 }
00097
00098 int Bernstein::nReferenceDOFsWithoutFacets(
00099 const CellType& maximalCellType,
00100 const CellType& cellType
00101 ) const
00102 {
00103 if (order_==0)
00104 {
00105 if (maximalCellType != cellType) return 0;
00106 return 1;
00107 }
00108
00109 switch(cellType)
00110 {
00111 case PointCell:
00112 return 1;
00113 case LineCell:
00114 return order_-1;
00115 case TriangleCell:
00116 if (order_ < 3) return 0;
00117 if (order_ >= 3) return (order_-1)*(order_-2)/2;
00118 break;
00119 case QuadCell:
00120 if (order_==1) return 0;
00121 if (order_==2) return 1;
00122 TEUCHOS_TEST_FOR_EXCEPTION(order_>2, std::runtime_error,
00123 "Bernstein order > 2 not implemented "
00124 "for quad cells");
00125 case TetCell:
00126 if (order_<=2) return 0;
00127 TEUCHOS_TEST_FOR_EXCEPTION(order_>2, std::runtime_error,
00128 "Bernstein order > 2 not implemented "
00129 "for tet cells");
00130 case BrickCell:
00131 if (order_<=1) return 0;
00132 if (order_==2) return 1;
00133 TEUCHOS_TEST_FOR_EXCEPTION(order_>2, std::runtime_error,
00134 "Bernstein order > 2 not implemented "
00135 "for brick cells");
00136 default:
00137 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00138 << cellType << " not implemented in Bernstein basis");
00139 return -1;
00140 }
00141 return -1;
00142 }
00143
00144 void Bernstein::getReferenceDOFs(
00145 const CellType& maximalCellType,
00146 const CellType& cellType,
00147 Array<Array<Array<int> > >& dofs) const
00148 {
00149 const int N=(order()+1)*(order()+2)/2;
00150 typedef Array<int> Aint;
00151 switch(cellType)
00152 {
00153 case PointCell:
00154 dofs.resize(1);
00155 dofs[0] = tuple<Aint>(tuple(0));
00156 return;
00157 case LineCell:
00158 dofs.resize(2);
00159 dofs[0] = tuple<Aint>(tuple(0), tuple(order()));
00160 dofs[1] = tuple<Aint>(makeRange(1, order()-1));
00161 return;
00162 case TriangleCell:
00163 {
00164 dofs.resize(3);
00165 dofs[0] = tuple<Aint>(tuple(0),tuple(N-order()-1),tuple(N-1));
00166 if (order()>1)
00167 {
00168
00169 Aint dof0(order()-1);
00170 Aint dof1;
00171 Aint dof2(order()-1);
00172
00173
00174 int inc = 2;
00175 dof0[0] = 1;
00176 for (int i=1;i<order()-1;i++)
00177 {
00178 dof0[i] = dof0[i-1]+inc;
00179 inc++;
00180 }
00181
00182
00183 dof1 = makeRange(N-order(),N-2);
00184
00185
00186 inc = 3;
00187 dof2[order()-2] = 2;
00188 for (int i=order()-3;i>=0;i--)
00189 {
00190 dof2[i] = dof2[i+1]+inc;
00191 inc++;
00192 }
00193
00194 dofs[1] = tuple<Aint>(dof1,dof2,dof0);
00195 }
00196 else
00197 {
00198 dofs[1] = tuple(Array<int>());
00199 }
00200 if (order()>2)
00201 {
00202 Array<int> internaldofs;
00203 int bfcur = 0;
00204
00205 for (int alpha1=order();alpha1>=0;alpha1--)
00206 {
00207 for (int alpha2=order()-alpha1;alpha2>=0;alpha2--)
00208 {
00209 int alpha3 = order()-alpha1-alpha2;
00210 if (alpha1*alpha2*alpha3>0) {
00211 internaldofs.append(bfcur);
00212 }
00213 bfcur++;
00214 }
00215 }
00216 dofs[2] = tuple(internaldofs);
00217 }
00218 else
00219 {
00220 dofs[2] = tuple(Array<int>());
00221 }
00222 return;
00223 }
00224 case TetCell:
00225 {
00226 }
00227 default:
00228 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00229 << cellType << " not implemented in Bernstein basis");
00230 }
00231 }
00232
00233
00234
00235 Array<int> Bernstein::makeRange(int low, int high)
00236 {
00237 if (high < low) return Array<int>();
00238
00239 Array<int> rtn(high-low+1);
00240 for (int i=0; i<rtn.length(); i++) rtn[i] = low+i;
00241 return rtn;
00242 }
00243
00244 void Bernstein::refEval(
00245 const CellType& cellType,
00246 const Array<Point>& pts,
00247 const SpatialDerivSpecifier& sds,
00248 Array<Array<Array<double> > >& result,
00249 int verbosity) const
00250 {
00251 TEUCHOS_TEST_FOR_EXCEPTION(!(sds.isPartial() || sds.isIdentity()),
00252 std::runtime_error,
00253 "cannot evaluate spatial derivative " << sds << " on Bernstein basis");
00254 const MultiIndex& deriv = sds.mi();
00255 typedef Array<double> Adouble;
00256 result.resize(1);
00257 result[0].resize(pts.length());
00258
00259 switch(cellType)
00260 {
00261 case PointCell:
00262 result[0] = tuple<Adouble>(tuple(1.0));
00263 return;
00264 case LineCell:
00265 for (int i=0; i<pts.length(); i++)
00266 {
00267 evalOnLine(pts[i], deriv, result[0][i]);
00268 }
00269 return;
00270 case TriangleCell:
00271 for (int i=0; i<pts.length(); i++)
00272 {
00273 evalOnTriangle(pts[i], deriv, result[0][i]);
00274 }
00275 return;
00276 case TetCell:
00277 for (int i=0; i<pts.length(); i++)
00278 {
00279 evalOnTet(pts[i], deriv, result[0][i]);
00280 }
00281 return;
00282 default:
00283 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00284 "Bernstein::refEval() unimplemented for cell type "
00285 << cellType);
00286
00287 }
00288 }
00289
00290
00291
00292 void Bernstein::evalOnLine(const Point& pt,
00293 const MultiIndex& deriv,
00294 Array<double>& result) const
00295 {
00296 ADReal x = ADReal(pt[0], 0, 1);
00297 ADReal one(1.0, 1);
00298
00299 result.resize(order()+1);
00300 Array<ADReal> tmp(result.length());
00301 Array<double> x0(order()+1);
00302
00303 if (order_ == 0)
00304 {
00305 tmp[0] = one;
00306 }
00307 else
00308 {
00309 double binom_cur = 1.0;
00310 for (int i=0; i<=order_; i++)
00311 {
00312 tmp[i] = one;
00313
00314 for (int j=0;j<order_-i;j++)
00315 {
00316 tmp[i] *= (1-x);
00317 }
00318 for (int j=0;j<i;j++)
00319 {
00320 tmp[i] *= x;
00321 }
00322 tmp[i] *= binom_cur;
00323 binom_cur *= double(order()-i) / double(i+1);
00324 }
00325 }
00326
00327 for (int i=0; i<tmp.length(); i++)
00328 {
00329 if (deriv.order()==0) result[i] = tmp[i].value();
00330 else result[i] = tmp[i].gradient()[0];
00331 }
00332 }
00333
00334 void Bernstein::evalOnTriangle(const Point& pt,
00335 const MultiIndex& deriv,
00336 Array<double>& result) const
00337
00338 {
00339 ADReal x = ADReal(pt[0], 0, 2);
00340 ADReal y = ADReal(pt[1], 1, 2);
00341 ADReal one(1.0, 2);
00342
00343 Array<ADReal> tmp;
00344
00345 SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
00346 << y.value());
00347
00348 if(order_==0) {
00349 result.resize(1);
00350 tmp.resize(1);
00351 tmp[0] = one;
00352 }
00353 else {
00354 int N = (order()+1)*(order()+2)/2;
00355 result.resize(N);
00356 tmp.resize(N);
00357
00358 ADReal b1 = 1.0 - x - y;
00359 ADReal b2 = x;
00360 ADReal b3 = y;
00361
00362
00363 int bfcur = 0;
00364
00365 for (int alpha1=order();alpha1>=0;alpha1--)
00366 {
00367 for (int alpha2 = order()-alpha1;alpha2>=0;alpha2--)
00368 {
00369 int alpha3 = order() - alpha1 - alpha2;
00370 tmp[bfcur] = one;
00371 for (int i=0;i<alpha1;i++)
00372 {
00373 tmp[bfcur] *= b1;
00374 }
00375 for (int i=0;i<alpha2;i++)
00376 {
00377 tmp[bfcur] *= b2;
00378 }
00379 for (int i=0;i<alpha3;i++)
00380 {
00381 tmp[bfcur] *= b3;
00382 }
00383 for (int i=2;i<=order();i++)
00384 {
00385 tmp[bfcur] *= (double) i;
00386 }
00387 for (int i=2;i<=alpha1;i++)
00388 {
00389 tmp[bfcur] /= (double) i;
00390 }
00391 for (int i=2;i<=alpha2;i++)
00392 {
00393 tmp[bfcur] /= (double) i;
00394 }
00395 for (int i=2;i<=alpha3;i++)
00396 {
00397 tmp[bfcur] /= (double) i;
00398 }
00399 bfcur++;
00400 }
00401 }
00402 }
00403
00404 for (int i=0; i<tmp.length(); i++)
00405 {
00406 if (deriv.order()==0) result[i] = tmp[i].value();
00407 else
00408 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00409 }
00410 }
00411
00412 void Bernstein::evalOnTet(const Point& pt,
00413 const MultiIndex& deriv,
00414 Array<double>& result) const
00415 {
00416 ADReal x = ADReal(pt[0], 0, 3);
00417 ADReal y = ADReal(pt[1], 1, 3);
00418 ADReal z = ADReal(pt[2], 2, 3);
00419 ADReal one(1.0, 3);
00420
00421 Array<ADReal> tmp(result.length());
00422
00423 if(order_==0)
00424 {
00425 tmp.resize(1);
00426 result.resize(1);
00427 tmp[0] = one;
00428 }
00429 else
00430 {
00431 }
00432
00433 for (int i=0; i<tmp.length(); i++)
00434 {
00435 if (deriv.order()==0) result[i] = tmp[i].value();
00436 else
00437 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00438 }
00439 }
00440
00441
00442
00443