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 "SundanceLegendre.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 Legendre::Legendre(int order):order_(order)
00044 {
00045
00046
00047
00048 if (order_ > 1) nrDOF_edge_ = order_ - 1;
00049 else nrDOF_edge_ = 0;
00050
00051 if (order_ > 3) nrDOF_face_ = ((order_-2)*(order_-3))/2;
00052 else nrDOF_face_ = 0;
00053
00054 nrDOF_brick_ = 0;
00055
00056
00057 }
00058
00059 bool Legendre::supportsCellTypePair(
00060 const CellType& maximalCellType,
00061 const CellType& cellType
00062 ) const
00063 {
00064 switch(maximalCellType)
00065 {
00066 case LineCell:
00067 switch(cellType)
00068 {
00069 case LineCell:
00070 case PointCell:
00071 return true;
00072 default:
00073 return false;
00074 }
00075 case QuadCell:
00076 switch(cellType)
00077 {
00078 case QuadCell:
00079 case LineCell:
00080 case PointCell:
00081 return true;
00082 default:
00083 return false;
00084 }
00085 case BrickCell:
00086 switch(cellType)
00087 {
00088 case BrickCell:
00089 case QuadCell:
00090 case LineCell:
00091 case PointCell:
00092 return true;
00093 default:
00094 return false;
00095 }
00096 default:
00097 return false;
00098 }
00099 }
00100
00101 void Legendre::print(std::ostream& os) const
00102 {
00103 os << "Legendre";
00104 }
00105
00106 int Legendre::nReferenceDOFsWithoutFacets(
00107 const CellType& maximalCellType,
00108 const CellType& cellType
00109 ) const
00110 {
00111 switch(cellType)
00112 {
00113 case PointCell:
00114 return 1;
00115 case LineCell:
00116 return nrDOF_edge_;
00117 case QuadCell:
00118 return nrDOF_face_;
00119 case BrickCell:
00120 return nrDOF_brick_;
00121 default:
00122 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "illegal combination of cell type and maximal cell type" );
00123 return -1;
00124 }
00125 }
00126
00127 Array<int> Legendre::makeRange(int low, int high)
00128 {
00129 if (high < low) return Array<int>();
00130
00131 Array<int> rtn(high-low+1);
00132 for (int i=0; i<rtn.length(); i++) rtn[i] = low+i;
00133 return rtn;
00134 }
00135
00136 void Legendre::getReferenceDOFs(
00137 const CellType& maximalCellType,
00138 const CellType& cellType,
00139 Array<Array<Array<int> > >& dofs) const
00140 {
00141
00142 typedef Array<int> Aint;
00143
00144
00145
00146 switch(cellType)
00147 {
00148 case PointCell:
00149 {
00150 dofs.resize(1);
00151 dofs[0] = tuple<Aint>(tuple(0));
00152 return;
00153 }
00154 break;
00155 case LineCell:
00156 {
00157 dofs.resize(2);
00158 dofs[0] = tuple<Aint>(tuple(0), tuple(1));
00159 if (nrDOF_edge_>0)
00160 {
00161 dofs[1] = tuple<Aint>(makeRange(2, 2+nrDOF_edge_-1));
00162 }
00163 else
00164 {
00165 dofs[1] = tuple(Array<int>());
00166 }
00167 return;
00168 }
00169 break;
00170 case QuadCell:
00171 {
00172 int offs = 0;
00173 dofs.resize(3);
00174
00175 dofs[0] = tuple<Aint>(tuple(0), tuple(1), tuple(2), tuple(3));
00176 offs = 4;
00177 if (nrDOF_edge_>0)
00178 {
00179 dofs[1].resize(4);
00180 dofs[1][0] = makeRange(offs, offs+nrDOF_edge_-1); offs += nrDOF_edge_;
00181 dofs[1][1] = makeRange(offs, offs+nrDOF_edge_-1); offs += nrDOF_edge_;
00182 dofs[1][2] = makeRange(offs, offs+nrDOF_edge_-1); offs += nrDOF_edge_;
00183 dofs[1][3] = makeRange(offs, offs+nrDOF_edge_-1); offs += nrDOF_edge_;
00184 }
00185 else
00186 {
00187 dofs[1] = tuple(Array<int>(), Array<int>(),
00188 Array<int>(), Array<int>());
00189 }
00190
00191 if (nrDOF_edge_>0)
00192 {
00193 dofs[2].resize(1);
00194 dofs[2][0] = makeRange(offs, offs+nrDOF_face_-1); offs += nrDOF_face_;
00195 }
00196 else
00197 {
00198 dofs[2] = tuple(Array<int>());
00199 }
00200
00201 }
00202 break;
00203 default:
00204 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00205 << cellType << " not implemented in Legendre basis");
00206 }
00207 }
00208
00209
00210 void Legendre::refEval(
00211 const CellType& cellType,
00212 const Array<Point>& pts,
00213 const SpatialDerivSpecifier& sds,
00214 Array<Array<Array<double> > >& result,
00215 int verbosity) const
00216 {
00217 TEUCHOS_TEST_FOR_EXCEPTION(!(sds.isPartial() || sds.isIdentity()),
00218 std::runtime_error,
00219 "cannot evaluate spatial derivative " << sds << " on Legendre basis");
00220 const MultiIndex& deriv = sds.mi();
00221 typedef Array<double> Adouble;
00222 result.resize(1);
00223 result[0].resize(pts.length());
00224
00225 switch(cellType)
00226 {
00227 case PointCell:
00228 result[0] = tuple<Adouble>(tuple(1.0));
00229 return;
00230 case LineCell:
00231 for (int i=0; i<pts.length(); i++)
00232 {
00233 evalOnLine(pts[i], deriv, result[0][i]);
00234 }
00235 return;
00236 case QuadCell:
00237 for (int i=0; i<pts.length(); i++)
00238 {
00239 evalOnQuad(pts[i], deriv, result[0][i]);
00240 }
00241 return;
00242 case BrickCell:
00243 for (int i=0; i<pts.length(); i++)
00244 {
00245 evalOnBrick(pts[i], deriv, result[0][i]);
00246 }
00247 return;
00248 default:
00249 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00250 "Legendre::refEval() unimplemented for cell type "
00251 << cellType);
00252
00253 }
00254 }
00255
00256
00257
00258 void Legendre::evalOnLine(const Point& pt,
00259 const MultiIndex& deriv,
00260 Array<double>& result) const
00261 {
00262 result.resize(2+nrDOF_edge_);
00263 ADReal x = ADReal(pt[0],0,1);
00264
00265 Array<ADReal> refAll(7);
00266
00267 refAll[0] = 1-x;
00268 refAll[1] = x;
00269 refAll[2] = 2.44948974278318 * ( (2*x-1)*(2*x-1) - 1 ) / 4;
00270 refAll[3] = 3.16227766016838 * ( (2*x-1)*(2*x-1) - 1 ) * (2*x-1) / 4;
00271 refAll[4] = 3.74165738677394 * ( 5*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 6*(2*x-1)*(2*x-1) + 1) / 16;
00272 refAll[5] = 4.24264068711929 * (2*x-1) * (7*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 10*(2*x-1)*(2*x-1) + 3) / 16;
00273 refAll[6] = 4.69041575982343 * (21*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) -
00274 35*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) + 15*(2*x-1)*(2*x-1) - 1) / 32;
00275
00276
00277 for (int i=0; i<result.length(); i++)
00278 {
00279 if (deriv.order()==0)
00280 {
00281 result[i] = refAll[i].value();
00282 }
00283 else
00284 {
00285 result[i] = refAll[i].gradient()[0];
00286 }
00287 }
00288
00289 return;
00290 }
00291
00292 void Legendre::evalOnQuad(const Point& pt,
00293 const MultiIndex& deriv,
00294 Array<double>& result) const
00295
00296 {
00297 result.resize( 4 + 4*nrDOF_edge_ + nrDOF_face_);
00298 ADReal x = ADReal(pt[0], 0, 2);
00299 ADReal y = ADReal(pt[1], 1, 2);
00300 ADReal one(1.0, 2);
00301
00302 Array<ADReal> refAllx(7);
00303 Array<ADReal> refAlly(7);
00304
00305 refAllx[0] = 1-x;
00306 refAllx[1] = x;
00307 refAllx[2] = 2.44948974278318 * ( (2*x-1)*(2*x-1) - 1 ) / 4;
00308 refAllx[3] = 3.16227766016838 * ( (2*x-1)*(2*x-1) - 1 ) * (2*x-1) / 4;
00309 refAllx[4] = 3.74165738677394 * ( 5*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 6*(2*x-1)*(2*x-1) + 1) / 16;
00310 refAllx[5] = 4.24264068711929 * (2*x-1) * (7*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) - 10*(2*x-1)*(2*x-1) + 3) / 16;
00311 refAllx[6] = 4.69041575982343 * (21*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) -
00312 35*(2*x-1)*(2*x-1)*(2*x-1)*(2*x-1) + 15*(2*x-1)*(2*x-1) - 1) / 32;
00313
00314 refAlly[0] = 1-y;
00315 refAlly[1] = y;
00316 refAlly[2] = 2.44948974278318 * ( (2*y-1)*(2*y-1) - 1 ) / 4;
00317 refAlly[3] = 3.16227766016838 * ( (2*y-1)*(2*y-1) - 1 ) * (2*y-1) / 4;
00318 refAlly[4] = 3.74165738677394 * ( 5*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) - 6*(2*y-1)*(2*y-1) + 1) / 16;
00319 refAlly[5] = 4.24264068711929 * (2*y-1) * (7*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) - 10*(2*y-1)*(2*y-1) + 3) / 16;
00320 refAlly[6] = 4.69041575982343 * (21*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) -
00321 35*(2*y-1)*(2*y-1)*(2*y-1)*(2*y-1) + 15*(2*y-1)*(2*y-1) - 1) / 32;
00322
00323 SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
00324 << y.value());
00325
00326 int sideIndex[4][2] = { {0,0} , {1,0} , {0,1} , {1,1}};
00327 int edgeIndex[4] = { 0 , 1 , 1 , 0};
00328 int edgeMultI[4] = { 0 , 0 , 1 , 1};
00329 int offs = 0;
00330 Array<ADReal> tmp(4 + 4*nrDOF_edge_ + nrDOF_face_);
00331
00332
00333 for (int i=0; i < 4 ; i++, offs++){
00334 tmp[offs] = refAllx[sideIndex[i][0]] * refAlly[sideIndex[i][1]];
00335 }
00336
00337
00338 for (int i=0; i < 4 ; i++){
00339
00340 if (edgeIndex[i] == 0){
00341 for (int j = 0 ; j < nrDOF_edge_ ; j++ , offs++){
00342 tmp[offs] = refAllx[2+j] * refAlly[edgeMultI[i]];
00343 }
00344 }
00345 else
00346 {
00347 for (int j = 0 ; j < nrDOF_edge_ ; j++ , offs++){
00348 tmp[offs] = refAllx[edgeMultI[i]] * refAlly[2+j];
00349 }
00350 }
00351 }
00352
00353
00354 if ( nrDOF_face_ > 0 ){
00355
00356 for (int hierarch = 0 ; hierarch < order_ - 3 ; hierarch++)
00357 {
00358
00359
00360 for (int i=0 ; i < hierarch+1 ; i++ , offs++)
00361 {
00362
00363 tmp[offs] = refAllx[2+i] * refAlly[2+(hierarch-i)];
00364 }
00365 }
00366 }
00367
00368
00369 for (int i=0; i<result.length(); i++)
00370 {
00371 if (deriv.order()==0) result[i] = tmp[i].value();
00372 else
00373 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00374 }
00375
00376 }
00377
00378 void Legendre::evalOnBrick(const Point& pt,
00379 const MultiIndex& deriv,
00380 Array<double>& result) const
00381 {
00382 ADReal x = ADReal(pt[0], 0, 3);
00383 ADReal y = ADReal(pt[1], 1, 3);
00384 ADReal z = ADReal(pt[2], 2, 3);
00385 ADReal one(1.0, 3);
00386
00387 TEUCHOS_TEST_FOR_EXCEPT(true);
00388 }