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 "SundanceElementIntegral.hpp"
00032 #include "SundanceOut.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "Teuchos_TimeMonitor.hpp"
00035
00036 using namespace Sundance;
00037 using namespace Teuchos;
00038
00039 using std::endl;
00040
00041 static Time& transCreationTimer()
00042 {
00043 static RCP<Time> rtn
00044 = TimeMonitor::getNewTimer("building integral transformation matrices");
00045 return *rtn;
00046 }
00047
00048
00049 bool& ElementIntegral::alwaysUseCofacets()
00050 {
00051 static bool rtn = true;
00052 return rtn;
00053 }
00054
00055 ElementIntegral::ElementIntegral(int spatialDim,
00056 const CellType& maxCellType,
00057 int dim,
00058 const CellType& cellType,
00059 bool isInternalBdry,
00060 const ParametrizedCurve& globalCurve,
00061 const Mesh& mesh,
00062 int verb)
00063 : setupVerb_(verb),
00064 integrationVerb_(0),
00065 transformVerb_(0),
00066 spatialDim_(spatialDim),
00067 dim_(dim),
00068 isInternalBdry_(isInternalBdry),
00069 nFacetCases_(1),
00070 testDerivOrder_(-1),
00071 nRefDerivTest_(-1),
00072 nNodesTest_(-1),
00073 unkDerivOrder_(-1),
00074 nRefDerivUnk_(-1),
00075 nNodesUnk_(-1),
00076 nNodes_(-1),
00077 order_(0),
00078 alpha_(),
00079 beta_(),
00080 cellType_(cellType),
00081 maxCellType_(maxCellType),
00082 evalCellType_(cellType),
00083 testBasis_(),
00084 unkBasis_(),
00085 globalCurve_(globalCurve),
00086 mesh_(mesh)
00087 {
00088 Tabs tab0;
00089 SUNDANCE_MSG2(setupVerb(), tab0 << "constructing 0-form ElementIntegral");
00090
00091
00092 if (alwaysUseCofacets() || (dim != spatialDim && !isInternalBdry))
00093 {
00094 evalCellType_ = maxCellType_;
00095 nFacetCases_ = numFacets(maxCellType, dim);
00096 }
00097 }
00098
00099 ElementIntegral::ElementIntegral(int spatialDim,
00100 const CellType& maxCellType,
00101 int dim,
00102 const CellType& cellType,
00103 const BasisFamily& testBasis,
00104 int alpha,
00105 int testDerivOrder,
00106 bool isInternalBdry,
00107 const ParametrizedCurve& globalCurve,
00108 const Mesh& mesh,
00109 int verb)
00110 : setupVerb_(verb),
00111 integrationVerb_(0),
00112 transformVerb_(0),
00113 spatialDim_(spatialDim),
00114 dim_(dim),
00115 isInternalBdry_(isInternalBdry),
00116 nFacetCases_(1),
00117 testDerivOrder_(testDerivOrder),
00118 nRefDerivTest_(ipow(spatialDim, testDerivOrder)),
00119 nNodesTest_(testBasis.nReferenceDOFsWithFacets(maxCellType, cellType)),
00120 unkDerivOrder_(-1),
00121 nRefDerivUnk_(-1),
00122 nNodesUnk_(-1),
00123 nNodes_(nNodesTest_),
00124 order_(1),
00125 alpha_(alpha),
00126 beta_(-1),
00127 cellType_(cellType),
00128 maxCellType_(maxCellType),
00129 evalCellType_(cellType),
00130 testBasis_(testBasis),
00131 unkBasis_(),
00132 globalCurve_(globalCurve),
00133 mesh_(mesh)
00134 {
00135 Tabs tab0(0);
00136 SUNDANCE_MSG2(setupVerb(), tab0 << "constructing 1-form ElementIntegral");
00137
00138
00139 bool okToRestrictTestToBdry = basisRestrictableToBoundary(testBasis);
00140
00141 Tabs tab1;
00142 SUNDANCE_MSG2(setupVerb(), tab1 << "dim=" << dim << " spatialDim=" << spatialDim);
00143 if (dim != spatialDim)
00144 {
00145 if (isInternalBdry)
00146 {
00147 TEUCHOS_TEST_FOR_EXCEPT(!okToRestrictTestToBdry);
00148 }
00149 if (alwaysUseCofacets() || testDerivOrder>0)
00150 {
00151 Tabs tab2;
00152 evalCellType_ = maxCellType_;
00153 nFacetCases_ = numFacets(maxCellType, dim);
00154 nNodesTest_ = testBasis.nReferenceDOFsWithFacets(maxCellType, maxCellType);
00155 SUNDANCE_MSG2(setupVerb(), tab2 << "nNodesTest=" << nNodesTest_);
00156 nNodes_ = nNodesTest_;
00157 TEUCHOS_TEST_FOR_EXCEPT(nNodes_ == 0);
00158 }
00159 else
00160 {
00161 TEUCHOS_TEST_FOR_EXCEPT(!okToRestrictTestToBdry);
00162 }
00163 }
00164
00165 SUNDANCE_MSG2(setupVerb(), tab1 << "nNodes=" << nNodes_);
00166 }
00167
00168
00169
00170 ElementIntegral::ElementIntegral(int spatialDim,
00171 const CellType& maxCellType,
00172 int dim,
00173 const CellType& cellType,
00174 const BasisFamily& testBasis,
00175 int alpha,
00176 int testDerivOrder,
00177 const BasisFamily& unkBasis,
00178 int beta,
00179 int unkDerivOrder,
00180 bool isInternalBdry,
00181 const ParametrizedCurve& globalCurve,
00182 const Mesh& mesh,
00183 int verb)
00184 : setupVerb_(verb),
00185 integrationVerb_(0),
00186 transformVerb_(0),
00187 spatialDim_(spatialDim),
00188 dim_(dim),
00189 isInternalBdry_(isInternalBdry),
00190 nFacetCases_(1),
00191 testDerivOrder_(testDerivOrder),
00192 nRefDerivTest_(ipow(spatialDim, testDerivOrder)),
00193 nNodesTest_(testBasis.nReferenceDOFsWithFacets(maxCellType, cellType)),
00194 unkDerivOrder_(unkDerivOrder),
00195 nRefDerivUnk_(ipow(spatialDim, unkDerivOrder)),
00196 nNodesUnk_(unkBasis.nReferenceDOFsWithFacets(maxCellType, cellType)),
00197 nNodes_(nNodesTest_*nNodesUnk_),
00198 order_(2),
00199 alpha_(alpha),
00200 beta_(beta),
00201 cellType_(cellType),
00202 maxCellType_(maxCellType),
00203 evalCellType_(cellType),
00204 testBasis_(testBasis),
00205 unkBasis_(unkBasis),
00206 globalCurve_(globalCurve),
00207 mesh_(mesh)
00208 {
00209 Tabs tab0(0);
00210 SUNDANCE_MSG2(setupVerb(), tab0
00211 << "***** constructing 2-form ElementIntegral");
00212
00213
00214 bool okToRestrictTestToBdry = basisRestrictableToBoundary(testBasis);
00215 bool okToRestrictUnkToBdry = basisRestrictableToBoundary(unkBasis);
00216
00217
00218 Tabs tab1;
00219 SUNDANCE_MSG2(setupVerb(), tab1 << "dim=" << dim << " spatialDim=" << spatialDim);
00220 if (dim != spatialDim)
00221 {
00222 if (isInternalBdry)
00223 {
00224 TEUCHOS_TEST_FOR_EXCEPT(!(okToRestrictTestToBdry && okToRestrictUnkToBdry));
00225 }
00226 if (alwaysUseCofacets() || testDerivOrder>0 || unkDerivOrder>0)
00227 {
00228 Tabs tab2;
00229 evalCellType_ = maxCellType_;
00230 nFacetCases_ = numFacets(maxCellType, dim);
00231 nNodesTest_ = testBasis.nReferenceDOFsWithFacets(maxCellType, maxCellType);
00232 SUNDANCE_MSG2(setupVerb(), tab2 << "nNodesTest=" << nNodesTest_);
00233 nNodesUnk_ = unkBasis.nReferenceDOFsWithFacets(maxCellType, maxCellType);
00234 SUNDANCE_MSG2(setupVerb(), tab2 << "nNodesUnk=" << nNodesUnk_);
00235 nNodes_ = nNodesTest_ * nNodesUnk_;
00236 TEUCHOS_TEST_FOR_EXCEPT(nNodes_ == 0);
00237 }
00238 else
00239 {
00240 TEUCHOS_TEST_FOR_EXCEPT(okToRestrictTestToBdry != okToRestrictUnkToBdry);
00241 }
00242 }
00243
00244 SUNDANCE_MSG2(setupVerb(), tab1 << "nNodes=" << nNodes_);
00245 }
00246
00247
00248 void ElementIntegral::setVerb(
00249 int integrationVerb,
00250 int transformVerb)
00251 {
00252 integrationVerb_ = integrationVerb;
00253 transformVerb_ = transformVerb;
00254 }
00255
00256 void ElementIntegral::describe(std::ostream& os) const
00257 {
00258 Tabs tab(0);
00259 bool hasTest = testDerivOrder() >= 0;
00260 bool hasUnk = unkDerivOrder() >= 0;
00261
00262 if (hasTest)
00263 {
00264 os << tab << "Test functions:" << std::endl;
00265 Tabs tab1;
00266 os << tab1 << "test basis = " << testBasis() << std::endl;
00267 os << tab1 << "test differentiation order = " << testDerivOrder() << std::endl;
00268 os << tab1 << "alpha = " << alpha() << std::endl;
00269 os << tab1 << "num test gradient components=" << nRefDerivTest() << std::endl;
00270 }
00271 if (hasUnk)
00272 {
00273 os << tab << "Unknown functions:" << std::endl;
00274 Tabs tab1;
00275 os << tab1 << "unk basis = " << unkBasis() << std::endl;
00276 os << tab1 << "unk differentiation order = " << unkDerivOrder() << std::endl;
00277 os << tab1 << "beta = " << beta() << std::endl;
00278 os << tab1 << "num unk gradient components=" << nRefDerivUnk() << std::endl;
00279 }
00280 os << tab << "Geometry:" << std::endl;
00281 Tabs tab1;
00282 os << tab1 << "cell type on which integral is defined: " << cellType()
00283 << std::endl;
00284 os << tab1 << "maximal cell type: " << maxCellType() << std::endl;
00285 os << tab1 << "cell type on which quad pts defined: " << evalCellType()
00286 << std::endl;
00287 os << tab << "Number of evaluation cases: " << nFacetCases() << std::endl;
00288 }
00289
00290
00291 void ElementIntegral::assertBilinearForm() const
00292 {
00293 TEUCHOS_TEST_FOR_EXCEPTION(testDerivOrder() < 0 || testDerivOrder() > 1,
00294 std::logic_error,
00295 "Test function derivative order=" << testDerivOrder()
00296 << " must be 0 or 1");
00297
00298 TEUCHOS_TEST_FOR_EXCEPTION(unkDerivOrder() < 0 || unkDerivOrder() > 1,
00299 std::logic_error,
00300 "Unknown function derivative order=" << unkDerivOrder()
00301 << " must be 0 or 1");
00302 }
00303
00304 void ElementIntegral::assertLinearForm() const
00305 {
00306 TEUCHOS_TEST_FOR_EXCEPTION(testDerivOrder() < 0 || testDerivOrder() > 1,
00307 std::logic_error,
00308 "Test function derivative order=" << testDerivOrder()
00309 << " must be 0 or 1");
00310 }
00311
00312
00313 void ElementIntegral::getQuad(const QuadratureFamily& quad,
00314 int evalCase, Array<Point>& quadPts, Array<double>& quadWeights) const
00315 {
00316 Tabs tab(0);
00317
00318 SUNDANCE_MSG2(setupVerb(), tab << "getting quad points for rule "
00319 << quad);
00320
00321 if (nFacetCases()==1)
00322 {
00323 quad.getPoints(cellType(), quadPts, quadWeights);
00324 }
00325 else
00326 {
00327 int dim = dimension(cellType());
00328 quad.getFacetPoints(maxCellType(), dim, evalCase,
00329 quadPts, quadWeights);
00330 }
00331
00332 if (setupVerb() >= 4)
00333 {
00334 Tabs tab1;
00335 SUNDANCE_MSG4(setupVerb(),
00336 tab1 << "quadrature points on ref cell are:");
00337 printQuad(Out::os(), quadPts, quadWeights);
00338 }
00339 }
00340
00341
00342
00343 Array<double>& ElementIntegral::G(int alpha)
00344 {
00345 static Array<Array<double> > rtn(3);
00346
00347 return rtn[alpha];
00348 }
00349
00350 Array<double>& ElementIntegral::G(int alpha, int beta)
00351 {
00352 static Array<Array<double> > rtn(9);
00353
00354 return rtn[3*alpha + beta];
00355 }
00356
00357 int& ElementIntegral::transformationMatrixIsValid(int alpha, int beta)
00358 {
00359 static Array<int> rtn(9, false);
00360 return rtn[3*alpha + beta];
00361 }
00362
00363 int& ElementIntegral::transformationMatrixIsValid(int alpha)
00364 {
00365 static Array<int> rtn(3, false);
00366 return rtn[alpha];
00367 }
00368
00369 void ElementIntegral::invalidateTransformationMatrices()
00370 {
00371 for (int i=0; i<3; i++)
00372 {
00373 transformationMatrixIsValid(i) = false;
00374 for (int j=0; j<3; j++)
00375 {
00376 transformationMatrixIsValid(i, j) = false;
00377 }
00378 }
00379 }
00380
00381
00382
00383
00384 int ElementIntegral::ipow(int base, int power)
00385 {
00386 int rtn = 1;
00387 for (int i=0; i<power; i++) rtn *= base;
00388 return rtn;
00389 }
00390
00391 void ElementIntegral
00392 ::createTwoFormTransformationMatrix(const CellJacobianBatch& JTrans,
00393 const CellJacobianBatch& JVol) const
00394 {
00395 TimeMonitor timer(transCreationTimer());
00396 Tabs tab;
00397
00398 int flops = 0;
00399
00400 int maxDim = JTrans.cellDim();
00401
00402
00403 if (testDerivOrder() == 1 && unkDerivOrder() == 1)
00404 {
00405 Tabs tab2;
00406 if (transformationMatrixIsValid(alpha(), beta())) return;
00407 transformationMatrixIsValid(alpha(), beta()) = true;
00408
00409 G(alpha(), beta()).resize(JTrans.numCells() * JTrans.cellDim() * JTrans.cellDim());
00410
00411 double* GPtr = &(G(alpha(),beta())[0]);
00412 int k = 0;
00413
00414 for (int c=0; c<JTrans.numCells(); c++)
00415 {
00416 static Array<double> invJ;
00417 JTrans.getInvJ(c, invJ);
00418 double detJ = fabs(JVol.detJ()[c]);
00419 for (int gamma=0; gamma<maxDim; gamma++)
00420 {
00421 for (int delta=0; delta<maxDim; delta++, k++)
00422 {
00423 GPtr[k] = detJ*invJ[alpha() + gamma*maxDim]
00424 * invJ[beta() + maxDim*delta];
00425 }
00426 }
00427 }
00428 flops = 2 * JTrans.numCells() * maxDim * maxDim + JTrans.numCells();
00429 }
00430
00431 else if (testDerivOrder() == 1 && unkDerivOrder() == 0)
00432 {
00433 if (transformationMatrixIsValid(alpha())) return;
00434 transformationMatrixIsValid(alpha()) = true;
00435
00436 G(alpha()).resize(JTrans.numCells() * JTrans.cellDim());
00437
00438 int k = 0;
00439 double* GPtr = &(G(alpha())[0]);
00440
00441 for (int c=0; c<JTrans.numCells(); c++)
00442 {
00443 static Array<double> invJ;
00444 JTrans.getInvJ(c, invJ);
00445 double detJ = fabs(JVol.detJ()[c]);
00446 for (int gamma=0; gamma<maxDim; gamma++,k++)
00447 {
00448 GPtr[k] = detJ*invJ[alpha() + maxDim * gamma];
00449 }
00450 }
00451 flops = JTrans.numCells() * maxDim + JTrans.numCells();
00452 }
00453
00454 else
00455 {
00456 if (transformationMatrixIsValid(beta())) return;
00457 transformationMatrixIsValid(beta()) = true;
00458
00459 G(beta()).resize(JTrans.numCells() * JTrans.cellDim());
00460
00461 int k = 0;
00462 double* GPtr = &(G(beta())[0]);
00463
00464 for (int c=0; c<JTrans.numCells(); c++)
00465 {
00466 static Array<double> invJ;
00467 JTrans.getInvJ(c, invJ);
00468 double detJ = fabs(JVol.detJ()[c]);
00469 for (int gamma=0; gamma<maxDim; gamma++,k++)
00470 {
00471 GPtr[k] = detJ*invJ[beta() + maxDim * gamma];
00472 }
00473 }
00474 flops = JTrans.numCells() * maxDim + JTrans.numCells();
00475 }
00476
00477 addFlops(flops);
00478 }
00479
00480
00481 void ElementIntegral
00482 ::createOneFormTransformationMatrix(const CellJacobianBatch& JTrans,
00483 const CellJacobianBatch& JVol) const
00484 {
00485 TimeMonitor timer(transCreationTimer());
00486 Tabs tab;
00487 SUNDANCE_MSG2(transformVerb(),
00488 tab << "ElementIntegral creating linear form trans matrices");
00489
00490 int maxDim = JTrans.cellDim();
00491
00492 if (transformationMatrixIsValid(alpha())) return;
00493 transformationMatrixIsValid(alpha()) = true;
00494
00495 int flops = JTrans.numCells() * maxDim + JTrans.numCells();
00496
00497 G(alpha()).resize(JTrans.numCells() * JTrans.cellDim());
00498
00499 int k = 0;
00500 double* GPtr = &(G(alpha())[0]);
00501
00502 for (int c=0; c<JTrans.numCells(); c++)
00503 {
00504 Array<double> invJ;
00505 JTrans.getInvJ(c, invJ);
00506 double detJ = fabs(JVol.detJ()[c]);
00507 for (int gamma=0; gamma<maxDim; gamma++, k++)
00508 {
00509 GPtr[k] = detJ*invJ[alpha() + maxDim * gamma];
00510 }
00511 }
00512
00513 addFlops(flops);
00514 }
00515