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 "SundanceBubble.hpp"
00032 #include "SundanceSpatialDerivSpecifier.hpp"
00033 #include "SundancePoint.hpp"
00034 #include "SundanceADReal.hpp"
00035 #include "PlayaExceptions.hpp"
00036 #include "SundanceObjectWithVerbosity.hpp"
00037 #include "SundanceOut.hpp"
00038
00039 using namespace Sundance;
00040 using namespace Sundance;
00041 using namespace Sundance;
00042 using namespace Teuchos;
00043
00044
00045
00046
00047 bool Bubble::supportsCellTypePair(
00048 const CellType& maximalCellType,
00049 const CellType& cellType
00050 ) const
00051 {
00052 switch(maximalCellType)
00053 {
00054 case LineCell:
00055 switch(cellType)
00056 {
00057 case LineCell:
00058 case PointCell:
00059 return true;
00060 default:
00061 return false;
00062 }
00063 case TriangleCell:
00064 switch(cellType)
00065 {
00066 case TriangleCell:
00067 case LineCell:
00068 case PointCell:
00069 return true;
00070 default:
00071 return false;
00072 }
00073 case TetCell:
00074 switch(cellType)
00075 {
00076 case TetCell:
00077 case TriangleCell:
00078 case LineCell:
00079 case PointCell:
00080 return true;
00081 default:
00082 return false;
00083 }
00084 default:
00085 return false;
00086 }
00087 }
00088
00089 void Bubble::print(std::ostream& os) const
00090 {
00091 os << "Bubble(" << order_ << ")";
00092 }
00093
00094 int Bubble::nReferenceDOFsWithoutFacets(
00095 const CellType& maximalCellType,
00096 const CellType& cellType
00097 ) const
00098 {
00099 if (maximalCellType==cellType) return 1;
00100 else return 0;
00101 }
00102
00103 void Bubble::getReferenceDOFs(
00104 const CellType& maxCellType,
00105 const CellType& cellType,
00106 Array<Array<Array<int> > >& dofs) const
00107 {
00108 typedef Array<int> Aint;
00109
00110 Aint z(1);
00111 z[0] = 0;
00112
00113 switch(cellType)
00114 {
00115 case PointCell:
00116 switch(maxCellType)
00117 {
00118 case LineCell:
00119 case TriangleCell:
00120 case TetCell:
00121 dofs.resize(1);
00122 dofs[0].resize(1);
00123 return;
00124 default:
00125 TEUCHOS_TEST_FOR_EXCEPT(1);
00126 }
00127 case LineCell:
00128 dofs.resize(2);
00129 dofs[0].resize(2);
00130 dofs[1].resize(1);
00131 if (maxCellType==LineCell)
00132 {
00133 dofs[1] = tuple<Aint>(z);
00134 }
00135 return;
00136 case TriangleCell:
00137 dofs.resize(3);
00138 dofs[0].resize(3);
00139 dofs[1].resize(3);
00140 dofs[2].resize(1);
00141 if (maxCellType==TriangleCell)
00142 {
00143 dofs[2]=tuple<Aint>(z);
00144 }
00145 return;
00146 case TetCell:
00147 dofs.resize(4);
00148 dofs[0].resize(4);
00149 dofs[1].resize(6);
00150 dofs[2].resize(4);
00151 dofs[3].resize(1);
00152 if (maxCellType==TetCell)
00153 {
00154 dofs[3]=tuple<Aint>(z);
00155 }
00156 return;
00157 default:
00158 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00159 << cellType << " not implemented in Bubble basis");
00160 }
00161 }
00162
00163
00164
00165
00166 void Bubble::refEval(
00167 const CellType& cellType,
00168 const Array<Point>& pts,
00169 const SpatialDerivSpecifier& sds,
00170 Array<Array<Array<double> > >& result,
00171 int verbosity) const
00172 {
00173 TEUCHOS_TEST_FOR_EXCEPTION(!(sds.isPartial() || sds.isIdentity()),
00174 std::runtime_error,
00175 "cannot evaluate spatial derivative " << sds << " on Bubble basis");
00176 const MultiIndex& deriv = sds.mi();
00177 typedef Array<double> Adouble;
00178 result.resize(1);
00179 result[0].resize(pts.length());
00180
00181 switch(cellType)
00182 {
00183 case PointCell:
00184 result[0] = tuple<Adouble>(tuple(1.0));
00185 return;
00186 case LineCell:
00187 for (int i=0; i<pts.length(); i++)
00188 {
00189 evalOnLine(pts[i], deriv, result[0][i]);
00190 }
00191 break;
00192 case TriangleCell:
00193 for (int i=0; i<pts.length(); i++)
00194 {
00195 evalOnTriangle(pts[i], deriv, result[0][i]);
00196 }
00197 break;
00198 case TetCell:
00199 for (int i=0; i<pts.length(); i++)
00200 {
00201 evalOnTet(pts[i], deriv, result[0][i]);
00202 }
00203 break;
00204 default:
00205 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00206 "Bubble::refEval() unimplemented for cell type "
00207 << cellType);
00208
00209 }
00210 }
00211
00212
00213
00214 void Bubble::evalOnLine(const Point& pt,
00215 const MultiIndex& deriv,
00216 Array<double>& result) const
00217 {
00218 ADReal x = ADReal(pt[0], 0, 1);
00219 ADReal one(1.0, 1);
00220
00221 result.resize(1);
00222 Array<ADReal> tmp(result.length());
00223
00224 int p = (order_+1)/2;
00225 if (order_==0) p=1;
00226 ADReal xp = one;
00227 for (int i=0; i<p; i++) xp = xp*x;
00228
00229 tmp[0] = xp*(1.0-xp);
00230
00231 for (int i=0; i<tmp.length(); i++)
00232 {
00233 if (deriv.order()==0) result[i] = tmp[i].value();
00234 else result[i] = tmp[i].gradient()[0];
00235 }
00236 }
00237
00238 void Bubble::evalOnTriangle(const Point& pt,
00239 const MultiIndex& deriv,
00240 Array<double>& result) const
00241
00242
00243
00244 {
00245 ADReal x = ADReal(pt[0], 0, 2);
00246 ADReal y = ADReal(pt[1], 1, 2);
00247 ADReal one(1.0, 2);
00248
00249 Array<ADReal> tmp(1);
00250 result.resize(1);
00251
00252
00253 SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
00254 << y.value());
00255
00256 int p = (order_+2)/3;
00257 if (order_==0) p=1;
00258
00259 ADReal xp = one;
00260 ADReal yp = one;
00261 for (int i=0; i<p; i++)
00262 {
00263 xp = xp*x;
00264 yp = yp*y;
00265 }
00266
00267 tmp[0] = xp*yp*(1.0-xp-yp)*std::pow(2.0,3*p);
00268
00269 for (int i=0; i<tmp.length(); i++)
00270 {
00271 SUNDANCE_OUT(this->verb() > 3,
00272 "tmp[" << i << "]=" << tmp[i].value()
00273 << " grad=" << tmp[i].gradient());
00274 if (deriv.order()==0) result[i] = tmp[i].value();
00275 else
00276 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00277 }
00278 }
00279
00280
00281 void Bubble::evalOnTet(const Point& pt,
00282 const MultiIndex& deriv,
00283 Array<double>& result) const
00284 {
00285 ADReal x = ADReal(pt[0], 0, 3);
00286 ADReal y = ADReal(pt[1], 1, 3);
00287 ADReal z = ADReal(pt[2], 2, 3);
00288 ADReal one(1.0, 3);
00289
00290 result.resize(1);
00291 Array<ADReal> tmp(result.length());
00292
00293 int p = (order_+3)/4;
00294 if (order_==0) p=1;
00295
00296 ADReal xp = one;
00297 ADReal yp = one;
00298 ADReal zp = one;
00299 for (int i=0; i<p; i++)
00300 {
00301 xp = xp*x;
00302 yp = yp*y;
00303 zp = zp*z;
00304 }
00305
00306 tmp[0] = (1.0-xp-yp-zp)*xp*yp*zp*std::pow(2.0,4*p);
00307
00308 for (int i=0; i<tmp.length(); i++)
00309 {
00310 if (deriv.order()==0) result[i] = tmp[i].value();
00311 else
00312 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00313 }
00314 }
00315