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
00032 #include "SundanceQuadratureFamily.hpp"
00033 #include "PlayaTabs.hpp"
00034
00035 using namespace Sundance;
00036 using namespace Teuchos;
00037 using Playa::Handle;
00038 using Playa::Handleable;
00039 using std::ios_base;
00040 using std::setw;
00041 using std::endl;
00042 using std::setprecision;
00043
00044 int QuadratureFamily::getNumPoints( const CellType & cellType ) const
00045 {
00046 const QuadratureFamilyBase* q
00047 = dynamic_cast<const QuadratureFamilyBase*>(ptr().get());
00048 return q->getNumPoints( cellType );
00049 }
00050
00051 void QuadratureFamily::getPoints(const CellType& cellType,
00052 Array<Point>& quadPoints,
00053 Array<double>& quadWeights) const
00054 {
00055 const QuadratureFamilyBase* q
00056 = dynamic_cast<const QuadratureFamilyBase*>(ptr().get());
00057
00058 TEUCHOS_TEST_FOR_EXCEPTION(q==0, std::logic_error,
00059 "QuadratureFamilyStub pointer" << toXML().toString()
00060 << " could not be cast to a QuadratureFamilyBase ptr");
00061
00062 q->getPoints(cellType, quadPoints, quadWeights);
00063 }
00064
00065 void QuadratureFamily::getAdaptedWeights(const CellType& cellType ,
00066 int cellDim,
00067 int celLID ,
00068 int facetIndex ,
00069 const Mesh& mesh ,
00070 const ParametrizedCurve& globalCurve ,
00071 Array<Point>& quadPoints ,
00072 Array<double>& quadWeights ,
00073 bool& isCut) const
00074 {
00075
00076 const QuadratureFamilyBase* q
00077 = dynamic_cast<const QuadratureFamilyBase*>(ptr().get());
00078
00079 TEUCHOS_TEST_FOR_EXCEPTION(q==0, std::logic_error,
00080 "QuadratureFamilyStub pointer" << toXML().toString()
00081 << " could not be cast to a QuadratureFamilyBase ptr");
00082
00083 q->getAdaptedWeights(cellType, cellDim , celLID , facetIndex ,
00084 mesh , globalCurve , quadPoints, quadWeights , isCut);
00085
00086 }
00087
00088 XMLObject QuadratureFamily::toXML() const
00089 {
00090 return ptr()->toXML();
00091 }
00092
00093 int QuadratureFamily::order() const
00094 {
00095 return ptr()->order();
00096 }
00097
00098
00099
00100 void QuadratureFamily::getFacetPoints(const CellType& cellType,
00101 int facetDim,
00102 int facetIndex,
00103 Array<Point>& quadPoints,
00104 Array<double>& quadWeights) const
00105 {
00106 switch(cellType)
00107 {
00108 case LineCell:
00109 getLineFacetQuad(facetDim, facetIndex, quadPoints, quadWeights);
00110 break;
00111 case TriangleCell:
00112 getTriangleFacetQuad(facetDim, facetIndex, quadPoints, quadWeights);
00113 break;
00114 case QuadCell:
00115 getQuadFacetQuad(facetDim, facetIndex, quadPoints, quadWeights);
00116 break;
00117 case TetCell:
00118 getTetFacetQuad(facetDim, facetIndex, quadPoints, quadWeights);
00119 break;
00120 case BrickCell:
00121 getBrickFacetQuad(facetDim, facetIndex, quadPoints, quadWeights);
00122 break;
00123 default:
00124 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00125 "getFacetPoints() not implemented for cell type "
00126 << cellType);
00127 }
00128 }
00129
00130
00131 void QuadratureFamily::getLineFacetQuad(int facetDim,
00132 int facetIndex,
00133 Array<Point>& quadPoints,
00134 Array<double>& quadWeights) const
00135 {
00136 TEUCHOS_TEST_FOR_EXCEPTION(facetDim > 0, std::runtime_error,
00137 "Invalid facet dimension " << facetDim
00138 << " in getLineFacetQuad()");
00139 TEUCHOS_TEST_FOR_EXCEPTION(facetIndex < 0 || facetIndex > 1, std::runtime_error,
00140 "Invalid facet index " << facetIndex
00141 << " in getLineFacetQuad()");
00142
00143 quadPoints.resize(1);
00144 quadWeights.resize(1);
00145 quadWeights[0] = 1.0;
00146 quadPoints[0] = Point((double) facetIndex);
00147 }
00148
00149
00150
00151
00152 void QuadratureFamily::getTriangleFacetQuad(int facetDim,
00153 int facetIndex,
00154 Array<Point>& quadPoints,
00155 Array<double>& quadWeights) const
00156 {
00157 TEUCHOS_TEST_FOR_EXCEPTION(facetDim > 1, std::runtime_error,
00158 "Invalid facet dimension " << facetDim
00159 << " in getTriangleFacetQuad()");
00160 TEUCHOS_TEST_FOR_EXCEPTION(facetIndex < 0 || facetIndex > 2, std::runtime_error,
00161 "Invalid facet index " << facetIndex
00162 << " in getTriangleFacetQuad()");
00163 if (facetDim==1)
00164 {
00165 Array<Point> facetPts;
00166 Array<double> facetWts;
00167 getPoints(LineCell, facetPts, facetWts);
00168 quadPoints.resize(facetPts.size());
00169 quadWeights.resize(facetWts.size());
00170 for (int i=0; i<facetPts.size(); i++)
00171 {
00172 if (facetIndex==2)
00173 {
00174 quadPoints[i] = Point(facetPts[i][0], 0.0);
00175 quadWeights[i] = facetWts[i];
00176 }
00177 else if (facetIndex==0)
00178 {
00179 quadPoints[i] = Point(1.0-facetPts[i][0], facetPts[i][0]);
00180 quadWeights[i] = facetWts[i];
00181 }
00182 else
00183 {
00184 quadPoints[i] = Point(0.0, facetPts[i][0]);
00185 quadWeights[i] = facetWts[i];
00186 }
00187 }
00188 }
00189 else
00190 {
00191 quadPoints.resize(1);
00192 quadWeights.resize(1);
00193 quadWeights[0] = 1.0;
00194 if (facetIndex==0) quadPoints[0] = Point(0.0, 0.0);
00195 else if (facetIndex==1) quadPoints[0] = Point(1.0, 0.0);
00196 else quadPoints[0] = Point(0.0, 1.0);
00197 }
00198 }
00199
00200 void QuadratureFamily::getQuadFacetQuad(int facetDim,
00201 int facetIndex,
00202 Array<Point>& quadPoints,
00203 Array<double>& quadWeights) const {
00204
00205 TEUCHOS_TEST_FOR_EXCEPTION(facetDim > 1, std::runtime_error,
00206 "Invalid facet dimension " << facetDim
00207 << " in getQuadFacetQuad()");
00208
00209 TEUCHOS_TEST_FOR_EXCEPTION(facetIndex < 0 || facetIndex > 3, std::runtime_error,
00210 "Invalid facet index " << facetIndex
00211 << " in getQuadFacetQuad()");
00212 if (facetDim==1)
00213 {
00214 Array<Point> facetPts;
00215 Array<double> facetWts;
00216 getPoints(LineCell, facetPts, facetWts);
00217 quadPoints.resize(facetPts.size());
00218 quadWeights.resize(facetWts.size());
00219 for (int i=0; i<facetPts.size(); i++)
00220 {
00221 if (facetIndex==0)
00222 {
00223 quadPoints[i] = Point(facetPts[i][0], 0.0);
00224 quadWeights[i] = facetWts[i];
00225 }
00226 else if (facetIndex==1)
00227 {
00228 quadPoints[i] = Point(0.0 , facetPts[i][0]);
00229 quadWeights[i] = facetWts[i];
00230 }
00231 else if (facetIndex==2)
00232 {
00233 quadPoints[i] = Point(1.0, facetPts[i][0]);
00234 quadWeights[i] = facetWts[i];
00235 }
00236 else
00237 {
00238 quadPoints[i] = Point(facetPts[i][0], 1.0);
00239 quadWeights[i] = facetWts[i];
00240 }
00241 }
00242 }
00243 else
00244 {
00245 quadPoints.resize(1);
00246 quadWeights.resize(1);
00247 quadWeights[0] = 1.0;
00248 if (facetIndex==0) quadPoints[0] = Point(0.0, 0.0);
00249 else if (facetIndex==1) quadPoints[0] = Point(1.0, 0.0);
00250 else if (facetIndex==2) quadPoints[0] = Point(0.0, 1.0);
00251 else quadPoints[0] = Point(1.0, 1.0);
00252 }
00253 }
00254
00255 void QuadratureFamily::getTetFacetQuad(int facetDim,
00256 int facetIndex,
00257 Array<Point>& quadPoints,
00258 Array<double>& quadWeights) const
00259 {
00260 TEUCHOS_TEST_FOR_EXCEPTION(facetDim > 2, std::runtime_error,
00261 "Invalid facet dimension " << facetDim
00262 << " in getTetFacetQuad()");
00263 TEUCHOS_TEST_FOR_EXCEPTION(facetIndex < 0 || facetIndex > 4, std::runtime_error,
00264 "Invalid facet index " << facetIndex
00265 << " in getTetFacetQuad()");
00266 if (facetDim==2)
00267 {
00268 Array<Point> facetPts;
00269 Array<double> facetWts;
00270 getPoints(TriangleCell, facetPts, facetWts);
00271 quadPoints.resize(facetPts.size());
00272 quadWeights.resize(facetWts.size());
00273 for (int i=0; i<facetPts.size(); i++)
00274 {
00275 double s = facetPts[i][0];
00276 double t = facetPts[i][1];
00277 double x,y,z;
00278 if (facetIndex==2)
00279 {
00280 x = s;
00281 y = 0.0;
00282 z = t;
00283 }
00284 else if (facetIndex==0)
00285 {
00286 x = 1.0-s-t;
00287 y = s;
00288 z = t;
00289 }
00290 else if (facetIndex==1)
00291 {
00292 x = 0.0;
00293 y = t;
00294 z = s;
00295 }
00296 else
00297 {
00298 x = s;
00299 y = 1.0-s-t;
00300 z = 0.0;
00301 }
00302 quadPoints[i] = Point(x, y, z);
00303 if (facetIndex != 1)
00304 quadWeights[i] = facetWts[i];
00305 else
00306 quadWeights[i] = facetWts[i] * sqrt(3.0);
00307
00308 }
00309 }
00310 else if (facetDim==0)
00311 {
00312 quadPoints.resize(1);
00313 quadWeights.resize(1);
00314 quadWeights[0] = 1.0;
00315 if (facetIndex==0) quadPoints[0] = Point(0.0, 0.0, 0.0);
00316 else if (facetIndex==1) quadPoints[0] = Point(1.0, 0.0, 0.0);
00317 else if (facetIndex==2) quadPoints[0] = Point(0.0, 1.0, 0.0);
00318 else quadPoints[0] = Point(0.0, 0.0, 1.0);
00319 }
00320 TEUCHOS_TEST_FOR_EXCEPT(facetDim==1);
00321 }
00322
00323
00324 void QuadratureFamily::getBrickFacetQuad(int facetDim,
00325 int facetIndex,
00326 Array<Point>& quadPoints,
00327 Array<double>& quadWeights) const
00328 {
00329 TEUCHOS_TEST_FOR_EXCEPTION(facetDim > 2, std::runtime_error,
00330 "Invalid facet dimension " << facetDim
00331 << " in getBrickFacetQuad()");
00332 if (facetDim==2)
00333 {
00334 TEUCHOS_TEST_FOR_EXCEPTION(facetIndex < 0 || facetIndex > 5, std::runtime_error,
00335 "Invalid facet index " << facetIndex
00336 << " in getBrickFacetQuad()");
00337 Array<Point> facetPts;
00338 Array<double> facetWts;
00339 getPoints(QuadCell, facetPts, facetWts);
00340 quadPoints.resize(facetPts.size());
00341 quadWeights.resize(facetWts.size());
00342 for (int i=0; i<facetPts.size(); i++)
00343 {
00344 if (facetIndex==0)
00345 {
00346 quadPoints[i] = Point(facetPts[i][0], facetPts[i][1] , 0.0);
00347 quadWeights[i] = facetWts[i];
00348 }
00349 else if (facetIndex==1)
00350 {
00351 quadPoints[i] = Point(facetPts[i][0], 0.0 , facetPts[i][1] );
00352 quadWeights[i] = facetWts[i];
00353 }
00354 else if (facetIndex==2)
00355 {
00356 quadPoints[i] = Point( 0.0 , facetPts[i][0] , facetPts[i][1] );
00357 quadWeights[i] = facetWts[i];
00358 }
00359 else if (facetIndex==3)
00360 {
00361 quadPoints[i] = Point( 1.0 , facetPts[i][0] , facetPts[i][1] );
00362 quadWeights[i] = facetWts[i];
00363 }
00364 else if (facetIndex==4)
00365 {
00366 quadPoints[i] = Point( facetPts[i][0] , 1.0 , facetPts[i][1] );
00367 quadWeights[i] = facetWts[i];
00368 }
00369 else
00370 {
00371 quadPoints[i] = Point( facetPts[i][0] , facetPts[i][1] , 1.0 );
00372 quadWeights[i] = facetWts[i];
00373 }
00374
00375 }
00376 }
00377 else if (facetDim==1)
00378 {
00379 TEUCHOS_TEST_FOR_EXCEPTION(facetIndex < 0 || facetIndex > 11, std::runtime_error,
00380 "Invalid facet index " << facetIndex
00381 << " in getBrickFacetQuad()");
00382 Array<Point> facetPts;
00383 Array<double> facetWts;
00384 getPoints(LineCell, facetPts, facetWts);
00385 quadPoints.resize(facetPts.size());
00386 quadWeights.resize(facetWts.size());
00387 for (int i=0; i<facetPts.size(); i++)
00388 {
00389 if (facetIndex==0)
00390 {
00391 quadPoints[i] = Point(facetPts[i][0], 0.0 , 0.0); quadWeights[i] = facetWts[i];
00392 }
00393 else if (facetIndex==1)
00394 { quadPoints[i] = Point( 0.0 , facetPts[i][0], 0.0 ); quadWeights[i] = facetWts[i]; }
00395 else if (facetIndex==2)
00396 { quadPoints[i] = Point( 0.0 , 0.0 , facetPts[i][0]); quadWeights[i] = facetWts[i]; }
00397 else if (facetIndex==3)
00398 { quadPoints[i] = Point( 1.0 , facetPts[i][0] , 0.0 ); quadWeights[i] = facetWts[i]; }
00399 else if (facetIndex==4)
00400 { quadPoints[i] = Point( 1.0 , 0.0 , facetPts[i][0] ); quadWeights[i] = facetWts[i]; }
00401 else if (facetIndex==5)
00402 { quadPoints[i] = Point( facetPts[i][0] , 1.0 , 0.0 ); quadWeights[i] = facetWts[i]; }
00403 else if (facetIndex==6)
00404 { quadPoints[i] = Point( 0.0 , 1.0 , facetPts[i][0]); quadWeights[i] = facetWts[i]; }
00405 else if (facetIndex==7)
00406 { quadPoints[i] = Point( 1.0 , 1.0 , facetPts[i][0]); quadWeights[i] = facetWts[i]; }
00407 else if (facetIndex==8)
00408 { quadPoints[i] = Point( facetPts[i][0] , 0.0 , 1.0); quadWeights[i] = facetWts[i]; }
00409 else if (facetIndex==9)
00410 { quadPoints[i] = Point( 0.0 , facetPts[i][0] , 1.0); quadWeights[i] = facetWts[i]; }
00411 else if (facetIndex==10)
00412 { quadPoints[i] = Point( 1.0 , facetPts[i][0] , 1.0); quadWeights[i] = facetWts[i]; }
00413 else
00414 { quadPoints[i] = Point( facetPts[i][0] , 1.0 , 1.0); quadWeights[i] = facetWts[i]; }
00415 }
00416 }
00417 else if (facetDim==0)
00418 {
00419 quadPoints.resize(1);
00420 quadWeights.resize(1);
00421 quadWeights[0] = 1.0;
00422 if (facetIndex==0) quadPoints[0] = Point(0.0, 0.0 , 0.0);
00423 else if (facetIndex==1) quadPoints[0] = Point(1.0, 0.0, 0.0);
00424 else if (facetIndex==2) quadPoints[0] = Point(0.0, 1.0, 0.0);
00425 else if (facetIndex==3) quadPoints[0] = Point(1.0, 1.0, 0.0);
00426 else if (facetIndex==4) quadPoints[0] = Point(0.0, 0.0, 1.0);
00427 else if (facetIndex==5) quadPoints[0] = Point(1.0, 0.0, 1.0);
00428 else if (facetIndex==6) quadPoints[0] = Point(0.0, 1.0, 1.0);
00429 else quadPoints[0] = Point(1.0, 1.0, 1.0);
00430 }
00431 }
00432
00433
00434 namespace Sundance
00435 {
00436 void printQuad(std::ostream& os,
00437 const Array<Point>& pts, const Array<double>& wgts)
00438 {
00439 Tabs tab(0);
00440 static Array<string> names = tuple<string>("x", "y", "z");
00441 TEUCHOS_TEST_FOR_EXCEPT(pts.size() != wgts.size());
00442
00443 TEUCHOS_TEST_FOR_EXCEPT(pts.size() < 1);
00444
00445 int dim = pts[0].dim();
00446
00447 int prec=10;
00448 int w = prec+5;
00449 ios_base::fmtflags oldFlags = os.flags();
00450 os.setf(ios_base::internal);
00451 os.setf(ios_base::showpoint);
00452
00453 os << tab << setw(5) << "i" << setw(w) << "w";;
00454
00455 for (int i=0; i<pts[0].dim(); i++)
00456 {
00457 os << setw(w) << names[i] ;
00458 }
00459 os << std::endl;
00460 os.setf(ios_base::right);
00461 for (int i=0; i<pts.size(); i++)
00462 {
00463 os << tab << setw(5) << i << setw(w) << setprecision(prec) << wgts[i];
00464 for (int d=0; d<dim; d++)
00465 {
00466 os << setw(w) << setprecision(prec) << pts[i][d];
00467 }
00468 os << std::endl;
00469 }
00470 os.flags(oldFlags);
00471 }
00472 }