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 "SundanceCubicHermite.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
00044 bool CubicHermite::supportsCellTypePair(
00045 const CellType& maximalCellType,
00046 const CellType& cellType
00047 ) const
00048 {
00049 switch(maximalCellType)
00050 {
00051 case LineCell:
00052 switch(cellType)
00053 {
00054 case LineCell:
00055 case PointCell:
00056 return true;
00057 default:
00058 return false;
00059 }
00060 case TriangleCell:
00061 switch(cellType)
00062 {
00063 case TriangleCell:
00064 case LineCell:
00065 case PointCell:
00066 return true;
00067 default:
00068 return false;
00069 }
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 default:
00082 return false;
00083 }
00084 }
00085
00086 void CubicHermite::print(std::ostream& os) const
00087 {
00088 os << "CubicHermite";
00089 }
00090
00091 int CubicHermite::nReferenceDOFsWithoutFacets(
00092 const CellType& maximalCellType,
00093 const CellType& cellType
00094 ) const
00095 {
00096 switch(maximalCellType)
00097 {
00098 case LineCell:
00099 switch(cellType)
00100 {
00101 case PointCell:
00102 return 2;
00103 case LineCell:
00104 return 0;
00105 default:
00106 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "illegal combination of cell type and maximal cell type" );
00107 return -1;
00108 }
00109 break;
00110 case TriangleCell:
00111 switch(cellType)
00112 {
00113 case PointCell:
00114 return 3;
00115 case LineCell:
00116 return 0;
00117 case TriangleCell:
00118 return 1;
00119 default:
00120 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "illegal combination of cell type and maximal cell type" );
00121 return -1;
00122 }
00123 break;
00124 default:
00125 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "illegal combination of cell type and maximal cell type" );
00126 return -1;
00127 }
00128
00129 }
00130
00131 void CubicHermite::getReferenceDOFs(
00132 const CellType& maximalCellType,
00133 const CellType& cellType,
00134 Array<Array<Array<int> > >& dofs) const
00135 {
00136 typedef Array<int> Aint;
00137 switch(cellType)
00138 {
00139 case PointCell:
00140 {
00141 dofs.resize(1);
00142 if (maximalCellType==LineCell)
00143 dofs[0] = tuple<Aint>(tuple(0,1));
00144 else if (maximalCellType==TriangleCell)
00145 dofs[0] = tuple<Aint>(tuple(0,1,2));
00146 else TEUCHOS_TEST_FOR_EXCEPT(1);
00147 return;
00148 }
00149 break;
00150 case LineCell:
00151 {
00152 dofs.resize(2);
00153 dofs[0].resize(2);
00154 if (maximalCellType==LineCell)
00155 {
00156 dofs[0][0].resize(2);
00157 dofs[0][0][0] = 0;
00158 dofs[0][0][1] = 1;
00159 dofs[0][1].resize(2);
00160 dofs[0][1][0] = 2;
00161 dofs[0][1][1] = 3;
00162 }
00163 else if (maximalCellType==TriangleCell)
00164 {
00165 dofs[0][0].resize(3);
00166 dofs[0][0][0] = 0;
00167 dofs[0][0][1] = 1;
00168 dofs[0][0][2] = 2;
00169 dofs[0][1].resize(3);
00170 dofs[0][1][0] = 3;
00171 dofs[0][1][1] = 4;
00172 dofs[0][1][1] = 5;
00173 }
00174 else
00175 {
00176 TEUCHOS_TEST_FOR_EXCEPT(1);
00177 }
00178 dofs[1].resize(1);
00179 dofs[1][0].resize(0);
00180 return;
00181 }
00182 break;
00183 case TriangleCell:
00184 {
00185 dofs.resize(3);
00186 dofs[0].resize(3);
00187 dofs[0][0] = tuple(0,1,2);
00188 dofs[0][1] = tuple(3,4,5);
00189 dofs[0][2] = tuple(6,7,8);
00190 dofs[1].resize(3);
00191 dofs[1][0].resize(0);
00192 dofs[1][1].resize(0);
00193 dofs[1][2].resize(0);
00194 dofs[2].resize(1);
00195 dofs[2][0].resize(1);
00196 dofs[2][0][0] = 9;
00197 }
00198 break;
00199 default:
00200 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00201 << cellType << " not implemented in CubicHermite basis");
00202 }
00203 }
00204
00205
00206 void CubicHermite::refEval(
00207 const CellType& cellType,
00208 const Array<Point>& pts,
00209 const SpatialDerivSpecifier& sds,
00210 Array<Array<Array<double> > >& result,
00211 int verbosity) const
00212 {
00213 TEUCHOS_TEST_FOR_EXCEPTION(!(sds.isPartial() || sds.isIdentity()),
00214 std::runtime_error,
00215 "cannot evaluate spatial derivative " << sds << " on CubicHermite basis");
00216 const MultiIndex& deriv = sds.mi();
00217 typedef Array<double> Adouble;
00218 result.resize(1);
00219 result[0].resize(pts.length());
00220
00221 switch(cellType)
00222 {
00223 case PointCell:
00224 result[0] = tuple<Adouble>(tuple(1.0));
00225 return;
00226 case LineCell:
00227 for (int i=0; i<pts.length(); i++)
00228 {
00229 evalOnLine(pts[i], deriv, result[0][i]);
00230 }
00231 return;
00232 case TriangleCell:
00233 for (int i=0; i<pts.length(); i++)
00234 {
00235 evalOnTriangle(pts[i], deriv, result[0][i]);
00236 }
00237 return;
00238 default:
00239 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00240 "CubicHermite::refEval() unimplemented for cell type "
00241 << cellType);
00242
00243 }
00244 }
00245
00246
00247
00248 void CubicHermite::evalOnLine(const Point& pt,
00249 const MultiIndex& deriv,
00250 Array<double>& result) const
00251 {
00252 result.resize(4);
00253 ADReal x = ADReal(pt[0],0,1);
00254 Array<ADReal> tmp(4);
00255
00256 tmp[0] = 1 + x * x * ( -3 + 2 * x );
00257 tmp[1] = x * ( 1 + (x - 2 ) * x );
00258 tmp[2] = ( 3 - 2*x ) * x * x;
00259 tmp[3] = (-1+x)*x*x;
00260
00261 for (int i=0; i<tmp.length(); i++)
00262 {
00263 if (deriv.order()==0)
00264 {
00265 result[i] = tmp[i].value();
00266 }
00267 else
00268 {
00269 result[i] = tmp[i].gradient()[0];
00270 }
00271 }
00272 return;
00273 }
00274
00275 void CubicHermite::evalOnTriangle(const Point& pt,
00276 const MultiIndex& deriv,
00277 Array<double>& result) const
00278
00279 {
00280 result.resize(10);
00281 ADReal x = ADReal(pt[0], 0, 2);
00282 ADReal y = ADReal(pt[1], 1, 2);
00283 ADReal one(1.0, 2);
00284
00285 Array<ADReal> tmp(10);
00286
00287 SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
00288 << y.value());
00289
00290 tmp[0] = 1 - 3*x*x + 2 * x*x*x - 13*x*y + 13*x*x*y - 3*y*y + 13 *x*y*y + 2 *y*y*y;
00291 tmp[1] = x - 2 *x*x + x*x*x - 3*x*y + 3*x*x*y + 2*x*y*y;
00292 tmp[2] = y - 3 *x *y + 2* x*x* y - 2* y*y + 3* x*y*y + y*y*y;
00293 tmp[3] = 3 * x*x - 2* x*x*x - 7* x* y + 7* x*x *y + 7* x*y*y;
00294 tmp[4] = -x*x + x*x*x + 2*x *y - 2* x*x* y - 2* x* y*y;
00295 tmp[5] = -x* y + 2* x*x* y + x* y*y;
00296 tmp[6] = -7* x* y + 7* x*x*y + 3* y*y + 7* x* y*y - 2* y*y*y;
00297 tmp[7] = -x *y + x*x* y + 2* x* y*y;
00298 tmp[8] = 2 *x *y - 2* x*x* y - y*y - 2* x* y*y + y*y*y;
00299 tmp[9] = 27* x *y - 27* x*x* y - 27* x* y*y;
00300
00301 for (int i=0; i<tmp.length(); i++)
00302 {
00303 if (deriv.order()==0) result[i] = tmp[i].value();
00304 else
00305 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00306 }
00307 }
00308
00309 void CubicHermite::evalOnTet(const Point& pt,
00310 const MultiIndex& deriv,
00311 Array<double>& result) const
00312 {
00313 ADReal x = ADReal(pt[0], 0, 3);
00314 ADReal y = ADReal(pt[1], 1, 3);
00315 ADReal z = ADReal(pt[2], 2, 3);
00316 ADReal one(1.0, 3);
00317
00318 TEUCHOS_TEST_FOR_EXCEPT(true);
00319 }
00320
00321 void CubicHermite::preApplyTransformation( const CellType &maxCellType ,
00322 const Mesh &mesh,
00323 const Array<int> &cellLIDs,
00324 const CellJacobianBatch& JVol,
00325 RCP<Array<double> >& A ) const
00326 {
00327 switch(maxCellType)
00328 {
00329 case TriangleCell:
00330 preApplyTransformationTriangle( mesh , cellLIDs, JVol , A );
00331 break;
00332 default:
00333 TEUCHOS_TEST_FOR_EXCEPT(1);
00334 break;
00335 }
00336 return;
00337 }
00338
00339 void CubicHermite::postApplyTransformation( const CellType &maxCellType ,
00340 const Mesh &mesh,
00341 const Array<int> &cellLIDs,
00342 const CellJacobianBatch& JVol,
00343 RCP<Array<double> >& A ) const
00344 {
00345 switch(maxCellType)
00346 {
00347 case TriangleCell:
00348 postApplyTransformationTriangle( mesh , cellLIDs, JVol , A );
00349 break;
00350 default:
00351 TEUCHOS_TEST_FOR_EXCEPT(1);
00352 break;
00353 }
00354 return;
00355 }
00356
00357 void CubicHermite::preApplyTransformationTranspose( const CellType &maxCellType ,
00358 const Mesh &mesh,
00359 const Array<int> &cellLIDs,
00360 const CellJacobianBatch& JVol ,
00361 Array<double> & A ) const
00362 {
00363 switch(maxCellType)
00364 {
00365 case TriangleCell:
00366 preApplyTransformationTransposeTriangle( mesh , cellLIDs, JVol , A );
00367 break;
00368 default:
00369 TEUCHOS_TEST_FOR_EXCEPT(1);
00370 break;
00371 }
00372 return;
00373 }
00374
00375
00376 void CubicHermite::preApplyTransformationTriangle( const Mesh &mesh,
00377 const Array<int> &cellLIDs,
00378 const CellJacobianBatch& JVol,
00379 RCP<Array<double> >& A) const
00380 {
00381 Array<double> &Anoptr = *A;
00382
00383 Array<double> cellVertH;
00384 getVertexH( mesh , cellLIDs , cellVertH );
00385
00386
00387
00388
00389
00390 const int numCols = Anoptr.size() / JVol.numCells() / 10;
00391
00392 for (int i=0;i<JVol.numCells();i++)
00393 {
00394 const int cell_start = i * numCols * 10;
00395 const double *invJ = JVol.jVals(i);
00396
00397 for (int j=0;j<numCols;j++)
00398 {
00399 const int col_start = cell_start + 10 * j;
00400 const double a1 = Anoptr[col_start + 1];
00401 const double a2 = Anoptr[col_start + 2];
00402 const double a4 = Anoptr[col_start + 4];
00403 const double a5 = Anoptr[col_start + 5];
00404 const double a7 = Anoptr[col_start + 7];
00405 const double a8 = Anoptr[col_start + 8];
00406 Anoptr[col_start+1] = (invJ[0]*a1 + invJ[1]*a2)/cellVertH[3*i];
00407 Anoptr[col_start+2] = (invJ[2]*a1 + invJ[3]*a2)/cellVertH[3*i];
00408 Anoptr[col_start+4] = (invJ[0]*a4 + invJ[1]*a5)/cellVertH[3*i+1];
00409 Anoptr[col_start+5] = (invJ[2]*a4 + invJ[3]*a5)/cellVertH[3*i+1];
00410 Anoptr[col_start+7] = (invJ[0]*a7 + invJ[1]*a8)/cellVertH[3*i+2];
00411 Anoptr[col_start+8] = (invJ[2]*a7 + invJ[3]*a8)/cellVertH[3*i+2];
00412 }
00413 }
00414
00415 }
00416
00417 void CubicHermite::postApplyTransformationTriangle( const Mesh &mesh,
00418 const Array<int> &cellLIDs ,
00419 const CellJacobianBatch& JVol,
00420 RCP<Array<double> >& A ) const
00421 {
00422 Array<double> &Anoptr = *A;
00423
00424 Array<double> cellVertH;
00425 getVertexH( mesh , cellLIDs , cellVertH );
00426
00427
00428 const int numRows = Anoptr.size() / 10 / JVol.numCells();
00429
00430 for (int i=0;i<JVol.numCells();i++)
00431 {
00432 const double *invJ = JVol.jVals(i);
00433
00434 const int cell_start = i * numRows * 10;
00435
00436 for (int j=0;j<numRows;j++)
00437 {
00438 const double a = Anoptr[cell_start + numRows + j];
00439 const double b = Anoptr[cell_start + 2*numRows + j];
00440 Anoptr[cell_start + numRows + j] = (invJ[0] * a + invJ[1] * b)/cellVertH[3*i];
00441 Anoptr[cell_start + 2*numRows + j] = (invJ[2] * a + invJ[3] * b)/cellVertH[3*i];
00442 }
00443
00444
00445 for (int j=0;j<numRows;j++)
00446 {
00447 const double a = Anoptr[cell_start + 4*numRows + j];
00448 const double b = Anoptr[cell_start + 5*numRows + j];
00449 Anoptr[cell_start + 4*numRows + j] = (invJ[0] * a + invJ[1] * b)/cellVertH[3*i+1];
00450 Anoptr[cell_start + 5*numRows + j] = (invJ[2] * a + invJ[3] * b)/cellVertH[3*i+1];
00451 }
00452
00453
00454 for (int j=0;j<numRows;j++)
00455 {
00456 const double a = Anoptr[cell_start + 7*numRows + j];
00457 const double b = Anoptr[cell_start + 8*numRows + j];
00458 Anoptr[cell_start + 7*numRows + j] = (invJ[0] * a + invJ[1] * b)/cellVertH[3*i+2];
00459 Anoptr[cell_start + 8*numRows + j] = (invJ[2] * a + invJ[3] * b)/cellVertH[3*i+2];
00460 }
00461
00462 }
00463
00464 return;
00465 }
00466
00467 void CubicHermite::preApplyTransformationTransposeTriangle( const Mesh &mesh,
00468 const Array<int> &cellLIDs,
00469 const CellJacobianBatch& JVol,
00470 Array<double>& A ) const
00471 {
00472
00473
00474
00475 const int numCols = A.size() / JVol.numCells() / 10;
00476
00477
00478 Array<double> cellVertH;
00479 getVertexH( mesh , cellLIDs , cellVertH );
00480
00481 for (int i=0;i<JVol.numCells();i++)
00482 {
00483 const int cell_start = i * numCols * 10;
00484 const double *invJ = JVol.jVals(i);
00485
00486 for (int j=0;j<numCols;j++)
00487 {
00488 const int col_start = cell_start + 10 * j;
00489 const double a1 = A[col_start + 1];
00490 const double a2 = A[col_start + 2];
00491 const double a4 = A[col_start + 4];
00492 const double a5 = A[col_start + 5];
00493 const double a7 = A[col_start + 7];
00494 const double a8 = A[col_start + 8];
00495 A[col_start+1] = (invJ[0]*a1 + invJ[2]*a2)/cellVertH[3*i];
00496 A[col_start+2] = (invJ[1]*a1 + invJ[3]*a2)/cellVertH[3*i];
00497 A[col_start+4] = (invJ[0]*a4 + invJ[2]*a5)/cellVertH[3*i+1];
00498 A[col_start+5] = (invJ[1]*a4 + invJ[3]*a5)/cellVertH[3*i+1];
00499 A[col_start+7] = (invJ[0]*a7 + invJ[2]*a8)/cellVertH[3*i+2];
00500 A[col_start+8] = (invJ[1]*a7 + invJ[3]*a8)/cellVertH[3*i+2];
00501 }
00502 }
00503
00504 }
00505
00506
00507
00508 void CubicHermite::getVertexH( const Mesh &mesh,
00509 const Array<int> &cellLIDs ,
00510 Array<double> &cellVertexH ) const
00511 {
00512 cellVertexH.resize(3*cellLIDs.size());
00513 const int N = mesh.numCells(2);
00514 const double h = 1.0 / sqrtl( N );
00515 for (int i=0;i<cellVertexH.size();i++)
00516 {
00517 cellVertexH[i] = h;
00518 }
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575 return;
00576
00577 }