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 "SundanceQuadratureIntegral.hpp"
00032 #include "SundanceGaussianQuadrature.hpp"
00033 #include "SundanceSpatialDerivSpecifier.hpp"
00034 #include "SundanceOut.hpp"
00035 #include "PlayaTabs.hpp"
00036 #include "Teuchos_TimeMonitor.hpp"
00037
00038 using namespace Sundance;
00039 using namespace Teuchos;
00040
00041 using std::endl;
00042 using std::setw;
00043 using std::setprecision;
00044
00045 extern "C"
00046 {
00047 int dgemm_(const char* transA, const char* transB,
00048 const int* M, const int *N, const int* K,
00049 const double* alpha,
00050 const double* A, const int* ldA,
00051 const double* B, const int* ldB,
00052 const double* beta,
00053 double* C, const int* ldC);
00054 }
00055
00056 static Time& quadrature0Timer()
00057 {
00058 static RCP<Time> rtn
00059 = TimeMonitor::getNewTimer("0-form quadrature");
00060 return *rtn;
00061 }
00062
00063 static Time& quadrature1Timer()
00064 {
00065 static RCP<Time> rtn
00066 = TimeMonitor::getNewTimer("1-form quadrature");
00067 return *rtn;
00068 }
00069
00070 static Time& quadrature2Timer()
00071 {
00072 static RCP<Time> rtn
00073 = TimeMonitor::getNewTimer("2-form quadrature");
00074 return *rtn;
00075 }
00076
00077
00078 QuadratureIntegral::QuadratureIntegral(int spatialDim,
00079 const CellType& maxCellType,
00080 int dim,
00081 const CellType& cellType,
00082 const QuadratureFamily& quad,
00083 bool isInternalBdry,
00084 const ParametrizedCurve& globalCurve,
00085 const Mesh& mesh,
00086 int verb)
00087 : QuadratureIntegralBase(spatialDim, maxCellType, dim, cellType, quad,
00088 isInternalBdry, globalCurve, mesh, verb),
00089 W_(),
00090 useSumFirstMethod_(true)
00091 {
00092 Tabs tab0(0);
00093
00094 SUNDANCE_MSG1(setupVerb(), tab0 << "QuadratureIntegral ctor for 0-form");
00095 if (setupVerb()) describe(Out::os());
00096
00097 SUNDANCE_MSG1(setupVerb(), tab0 << "quadrature family=" << quad);
00098
00099
00100 W_.resize(nFacetCases());
00101 quadPts_.resize(nFacetCases());
00102 quadWeights_.resize(nFacetCases());
00103 quad_ = quad;
00104
00105 for (int fc=0; fc<nFacetCases(); fc++)
00106 {
00107 Tabs tab2;
00108 SUNDANCE_MSG1(setupVerb(), tab2 << "facet case=" << fc
00109 << " of " << nFacetCases());
00110 Tabs tab3;
00111
00112 Array<double> quadWeights;
00113 Array<Point> quadPts;
00114 getQuad(quad, fc, quadPts, quadWeights);
00115 nQuad_ = quadPts.size();
00116 quadPts_[fc] = quadPts;
00117 quadWeights_[fc] = quadWeights;
00118
00119 W_[fc].resize(nQuad());
00120
00121 for (int q=0; q<nQuad(); q++)
00122 {
00123 W_[fc][q] = quadWeights[q];
00124 }
00125 }
00126 }
00127
00128
00129 QuadratureIntegral::QuadratureIntegral(int spatialDim,
00130 const CellType& maxCellType,
00131 int dim,
00132 const CellType& cellType,
00133 const BasisFamily& testBasis,
00134 int alpha,
00135 int testDerivOrder,
00136 const QuadratureFamily& quad,
00137 bool isInternalBdry,
00138 const ParametrizedCurve& globalCurve,
00139 const Mesh& mesh,
00140 int verb)
00141 : QuadratureIntegralBase(spatialDim, maxCellType, dim, cellType,
00142 testBasis, alpha, testDerivOrder, quad , isInternalBdry, globalCurve , mesh, verb),
00143 W_(),
00144 useSumFirstMethod_(true)
00145 {
00146 Tabs tab0;
00147
00148 SUNDANCE_MSG1(setupVerb(), tab0 << "QuadratureIntegral ctor for 1-form");
00149 if (setupVerb()) describe(Out::os());
00150 assertLinearForm();
00151
00152
00153 SUNDANCE_MSG1(setupVerb(), tab0 << "quadrature family=" << quad);
00154
00155 quad_ = quad;
00156 W_.resize(nFacetCases());
00157
00158 quadPts_.resize(nFacetCases());
00159 quadWeights_.resize(nFacetCases());
00160 W_ACI_F1_.resize(nFacetCases());
00161
00162 for (int fc=0; fc<nFacetCases(); fc++)
00163 {
00164 Tabs tab2;
00165
00166 SUNDANCE_MSG1(setupVerb(), tab2 << "facet case=" << fc
00167 << " of " << nFacetCases());
00168
00169
00170 Array<double> quadWeights;
00171 Array<Point> quadPts;
00172 getQuad(quad, fc, quadPts, quadWeights);
00173 nQuad_ = quadPts.size();
00174 quadPts_[fc] = quadPts;
00175 quadWeights_[fc] = quadWeights;
00176
00177 W_[fc].resize(nQuad() * nRefDerivTest() * nNodesTest());
00178
00179 SUNDANCE_MSG1(setupVerb(), tab2 << "num nodes for test function " << nNodesTest());
00180
00181 Array<Array<Array<Array<double> > > > testBasisVals(nRefDerivTest());
00182
00183 for (int r=0; r<nRefDerivTest(); r++)
00184 {
00185 Tabs tab3;
00186 SUNDANCE_MSG1(setupVerb(), tab3
00187 << "evaluating basis functions for ref deriv direction " << r);
00188 MultiIndex mi;
00189 testBasisVals[r].resize(testBasis.dim());
00190 if (testDerivOrder==1) mi[r] = 1;
00191 SpatialDerivSpecifier deriv(mi);
00192 testBasis.refEval(evalCellType(), quadPts, deriv,
00193 testBasisVals[r], setupVerb());
00194 }
00195
00196
00197
00198 int vecComp = 0;
00199 W_ACI_F1_[fc].resize(nQuad());
00200 for (int q=0; q<nQuad(); q++)
00201 {
00202 W_ACI_F1_[fc][q].resize(nRefDerivTest());
00203 for (int t=0; t<nRefDerivTest(); t++)
00204 {
00205 W_ACI_F1_[fc][q][t].resize(nNodesTest());
00206 for (int nt=0; nt<nNodesTest(); nt++)
00207 {
00208 wValue(fc, q, t, nt)
00209 = chop(quadWeights[q] * testBasisVals[t][vecComp][q][nt]) ;
00210 W_ACI_F1_[fc][q][t][nt] = chop(testBasisVals[t][vecComp][q][nt]);
00211 }
00212 }
00213 }
00214
00215 addFlops(2*nQuad()*nRefDerivTest()*nNodesTest());
00216 }
00217 }
00218
00219
00220
00221
00222 QuadratureIntegral::QuadratureIntegral(int spatialDim,
00223 const CellType& maxCellType,
00224 int dim,
00225 const CellType& cellType,
00226 const BasisFamily& testBasis,
00227 int alpha,
00228 int testDerivOrder,
00229 const BasisFamily& unkBasis,
00230 int beta,
00231 int unkDerivOrder,
00232 const QuadratureFamily& quad,
00233 bool isInternalBdry,
00234 const ParametrizedCurve& globalCurve,
00235 const Mesh& mesh,
00236 int verb)
00237 : QuadratureIntegralBase(spatialDim, maxCellType, dim, cellType,
00238 testBasis, alpha, testDerivOrder,
00239 unkBasis, beta, unkDerivOrder, quad , isInternalBdry, globalCurve , mesh , verb),
00240 W_(),
00241 useSumFirstMethod_(true)
00242 {
00243 Tabs tab0;
00244
00245 SUNDANCE_MSG1(setupVerb(), tab0 << "QuadratureIntegral ctor for 2-form");
00246 if (setupVerb()) describe(Out::os());
00247 assertBilinearForm();
00248 quad_ = quad;
00249
00250 W_.resize(nFacetCases());
00251 W_ACI_F2_.resize(nFacetCases());
00252
00253 quadPts_.resize(nFacetCases());
00254 quadWeights_.resize(nFacetCases());
00255
00256 for (int fc=0; fc<nFacetCases(); fc++)
00257 {
00258
00259 Array<double> quadWeights;
00260 Array<Point> quadPts;
00261 getQuad(quad, fc, quadPts, quadWeights);
00262
00263 quadPts_[fc] = quadPts;
00264 quadWeights_[fc] = quadWeights;
00265 nQuad_ = quadPts.size();
00266
00267 W_[fc].resize(nQuad() * nRefDerivTest() * nNodesTest()
00268 * nRefDerivUnk() * nNodesUnk());
00269
00270
00271
00272 Array<Array<Array<Array<double> > > > testBasisVals(nRefDerivTest());
00273 Array<Array<Array<Array<double> > > > unkBasisVals(nRefDerivUnk());
00274
00275
00276 for (int r=0; r<nRefDerivTest(); r++)
00277 {
00278 testBasisVals[r].resize(testBasis.dim());
00279 MultiIndex mi;
00280 if (testDerivOrder==1) mi[r] = 1;
00281 SpatialDerivSpecifier deriv(mi);
00282 testBasis.refEval(evalCellType(), quadPts, deriv,
00283 testBasisVals[r], setupVerb());
00284 }
00285
00286 for (int r=0; r<nRefDerivUnk(); r++)
00287 {
00288 unkBasisVals[r].resize(unkBasis.dim());
00289 MultiIndex mi;
00290 if (unkDerivOrder==1) mi[r] = 1;
00291 SpatialDerivSpecifier deriv(mi);
00292 unkBasis.refEval(evalCellType(),
00293 quadPts, deriv, unkBasisVals[r], setupVerb());
00294 }
00295
00296
00297 int vecComp = 0;
00298
00299 W_ACI_F2_[fc].resize(nQuad());
00300 for (int q=0; q<nQuad(); q++)
00301 {
00302 W_ACI_F2_[fc][q].resize(nRefDerivTest());
00303 for (int t=0; t<nRefDerivTest(); t++)
00304 {
00305 W_ACI_F2_[fc][q][t].resize(nNodesTest());
00306 for (int nt=0; nt<nNodesTest(); nt++)
00307 {
00308 W_ACI_F2_[fc][q][t][nt].resize(nRefDerivUnk());
00309 for (int u=0; u<nRefDerivUnk(); u++)
00310 {
00311 W_ACI_F2_[fc][q][t][nt][u].resize(nNodesUnk());
00312 for (int nu=0; nu<nNodesUnk(); nu++)
00313 {
00314 wValue(fc, q, t, nt, u, nu)
00315 = chop(quadWeights[q] * testBasisVals[t][vecComp][q][nt]
00316 * unkBasisVals[u][vecComp][q][nu]);
00317 W_ACI_F2_[fc][q][t][nt][u][nu] =
00318 chop(testBasisVals[t][vecComp][q][nt] * unkBasisVals[u][vecComp][q][nu]);
00319 }
00320 }
00321 }
00322 }
00323 }
00324
00325 addFlops(3*nQuad()*nRefDerivTest()*nNodesTest()*nRefDerivUnk()*nNodesUnk()
00326 + W_[fc].size());
00327 for (int i=0; i<W_[fc].size(); i++) W_[fc][i] = chop(W_[fc][i]);
00328
00329 }
00330
00331 }
00332
00333
00334 void QuadratureIntegral::transformZeroForm(const CellJacobianBatch& JTrans,
00335 const CellJacobianBatch& JVol,
00336 const Array<int>& isLocalFlag,
00337 const Array<int>& facetIndex,
00338 const RCP<Array<int> >& cellLIDs,
00339 const double* const coeff,
00340 RCP<Array<double> >& A) const
00341 {
00342 TimeMonitor timer(quadrature0Timer());
00343 Tabs tabs;
00344 TEUCHOS_TEST_FOR_EXCEPTION(order() != 0, std::logic_error,
00345 "QuadratureIntegral::transformZeroForm() called "
00346 "for form of order " << order());
00347
00348 TEUCHOS_TEST_FOR_EXCEPTION( (int) isLocalFlag.size() != 0
00349 && (int) isLocalFlag.size() != JVol.numCells(),
00350 std::runtime_error,
00351 "mismatch between isLocalFlag.size()=" << isLocalFlag.size()
00352 << " and JVol.numCells()=" << JVol.numCells());
00353
00354 bool checkLocalFlag = (int) isLocalFlag.size() != 0;
00355
00356 const Array<int>* cellLID = cellLIDs.get();
00357
00358 SUNDANCE_MSG1(integrationVerb(), tabs << "doing zero form by quadrature");
00359 double& a = (*A)[0];
00360 SUNDANCE_MSG5(integrationVerb(), tabs << "input A=");
00361 if (integrationVerb() >= 5) writeTable(Out::os(), tabs, *A, 6);
00362 double* coeffPtr = (double*) coeff;
00363
00364 int flops = 0;
00365 if (nFacetCases()==1)
00366 {
00367 if (globalCurve().isCurveValid())
00368 {
00369
00370 Array<double> quadWeightsTmp = quadWeights_[0];
00371 Array<Point> quadPointsTmp = quadPts_[0];
00372 bool isCutByCurve;
00373
00374 const Array<double>& w = W_[0];
00375 for (int c=0; c<JVol.numCells(); c++)
00376 {
00377 if (checkLocalFlag && !isLocalFlag[c])
00378 {
00379 coeffPtr += nQuad();
00380 continue;
00381 }
00382 double detJ = fabs(JVol.detJ()[c]);
00383 int fc = 0;
00384
00385 quadWeightsTmp = quadWeights_[fc];
00386 quadPointsTmp = quadPts_[fc];
00387 quad_.getAdaptedWeights(cellType(), dim(), (*cellLID)[c], fc ,mesh(),
00388 globalCurve(), quadPointsTmp, quadWeightsTmp, isCutByCurve);
00389
00390 if (isCutByCurve)
00391 {
00392 for (int q=0; q<nQuad(); q++, coeffPtr++)
00393 {
00394 a += quadWeightsTmp[q]*(*coeffPtr)*detJ;
00395 }
00396 flops += 3*nQuad();
00397 }
00398 else
00399 {
00400 for (int q=0; q<nQuad(); q++, coeffPtr++)
00401 {
00402 a += w[q]*(*coeffPtr)*detJ;
00403 }
00404 flops += 3*nQuad();
00405 }
00406
00407 }
00408 }
00409 else
00410 {
00411 const Array<double>& w = W_[0];
00412 for (int c=0; c<JVol.numCells(); c++)
00413 {
00414 if (checkLocalFlag && !isLocalFlag[c])
00415 {
00416 coeffPtr += nQuad();
00417 continue;
00418 }
00419 double detJ = fabs(JVol.detJ()[c]);
00420
00421 for (int q=0; q<nQuad(); q++, coeffPtr++)
00422 {
00423 a += w[q]*(*coeffPtr)*detJ;
00424 }
00425 flops += 3*nQuad();
00426 }
00427 }
00428 }
00429 else
00430 {
00431 if (globalCurve().isCurveValid())
00432 {
00433 Array<double> quadWeightsTmp = quadWeights_[0];
00434 Array<Point> quadPointsTmp = quadPts_[0];
00435 bool isCutByCurve;
00436
00437 for (int c=0; c<JVol.numCells(); c++)
00438 {
00439 if (checkLocalFlag && !isLocalFlag[c])
00440 {
00441 coeffPtr += nQuad();
00442 continue;
00443 }
00444 double detJ = fabs(JVol.detJ()[c]);
00445 int fc = facetIndex[c];
00446 const Array<double>& w = W_[facetIndex[c]];
00447
00448
00449 quadWeightsTmp = quadWeights_[fc];
00450 quadPointsTmp = quadPts_[fc];
00451 quad_.getAdaptedWeights(cellType(), dim(), (*cellLID)[c], fc ,mesh(),
00452 globalCurve(), quadPointsTmp, quadWeightsTmp, isCutByCurve);
00453
00454 if (isCutByCurve){
00455 for (int q=0; q<nQuad(); q++, coeffPtr++)
00456 a += quadWeightsTmp[q]*(*coeffPtr)*detJ;
00457 flops += 3*nQuad();
00458 }
00459 else{
00460 for (int q=0; q<nQuad(); q++, coeffPtr++)
00461 a += w[q]*(*coeffPtr)*detJ;
00462 flops += 3*nQuad();
00463 }
00464 }
00465 }
00466 else
00467 {
00468 for (int c=0; c<JVol.numCells(); c++)
00469 {
00470 Tabs tab1;
00471 SUNDANCE_MSG4(integrationVerb(), tab1 << "cell #" << c
00472 << " of " << JVol.numCells());
00473 if (checkLocalFlag && !isLocalFlag[c])
00474 {
00475 Tabs tab2;
00476 SUNDANCE_MSG4(integrationVerb(), tab2
00477 << "skipped -- is not local cell");
00478 coeffPtr += nQuad();
00479 continue;
00480 }
00481 double detJ = fabs(JVol.detJ()[c]);
00482 const Array<double>& w = W_[facetIndex[c]];
00483
00484 for (int q=0; q<nQuad(); q++, coeffPtr++)
00485 {
00486 Tabs tab3;
00487 SUNDANCE_MSG4(integrationVerb(),
00488 tab3 << "q=" << q << "\t w=" << w[q]
00489 << "\t\t coeff=" << *coeffPtr << "\t\t |J|=" << detJ);
00490 a += w[q]*(*coeffPtr)*detJ;
00491 }
00492 flops += 3*nQuad();
00493 }
00494 }
00495 }
00496 SUNDANCE_MSG5(integrationVerb(), tabs << "output A = ");
00497 if (integrationVerb() >= 5) writeTable(Out::os(), tabs, *A, 6);
00498
00499 addFlops(flops);
00500 }
00501
00502
00503 void QuadratureIntegral::transformOneForm(const CellJacobianBatch& JTrans,
00504 const CellJacobianBatch& JVol,
00505 const Array<int>& facetIndex,
00506 const RCP<Array<int> >& cellLIDs,
00507 const double* const coeff,
00508 RCP<Array<double> >& A) const
00509 {
00510 TimeMonitor timer(quadrature1Timer());
00511 Tabs tabs;
00512 TEUCHOS_TEST_FOR_EXCEPTION(order() != 1, std::logic_error,
00513 "QuadratureIntegral::transformOneForm() called for form "
00514 "of order " << order());
00515 SUNDANCE_MSG2(integrationVerb(), tabs << "doing one form by quadrature");
00516 int flops = 0;
00517
00518 const Array<int>* cellLID = cellLIDs.get();
00519
00520
00521
00522
00523 if (testDerivOrder() == 0)
00524 {
00525 double* aPtr = &((*A)[0]);
00526 SUNDANCE_MSG5(integrationVerb(), tabs << "input A = ");
00527 if (integrationVerb() >= 5) writeTable(Out::os(), tabs, *A, 6);
00528
00529 double* coeffPtr = (double*) coeff;
00530 int offset = 0 ;
00531
00532 if (globalCurve().isCurveValid())
00533 {
00534
00535 Array<double> quadWeightsTmp = quadWeights_[0];
00536 Array<Point> quadPointsTmp = quadPts_[0];
00537 bool isCutByCurve;
00538
00539 for (int c=0; c<JVol.numCells(); c++, offset+=nNodes())
00540 {
00541 Tabs tab2;
00542 double detJ = fabs(JVol.detJ()[c]);
00543 int fc = 0;
00544 if (nFacetCases() != 1) fc = facetIndex[c];
00545 const Array<double>& w = W_[fc];
00546 SUNDANCE_MSG4(integrationVerb(), tab2 << "c=" << c << " detJ=" << detJ);
00547
00548 quadWeightsTmp = quadWeights_[fc];
00549 quadPointsTmp = quadPts_[fc];
00550
00551 quad_.getAdaptedWeights(cellType(), dim(), (*cellLID)[c] , fc ,
00552 mesh() , globalCurve() , quadPointsTmp , quadWeightsTmp , isCutByCurve );
00553 if (isCutByCurve){
00554 Array<double> wi;
00555 wi.resize(nQuad() * nNodes());
00556 for (int ii = 0 ; ii < wi.size() ; ii++ ) { wi[ii] = 0.0; }
00557 for (int n = 0 ; n < nNodes() ; n++)
00558 for (int q=0 ; q < quadWeightsTmp.size() ; q++)
00559
00560 wi[n + nNodes()*q] +=
00561 chop(quadWeightsTmp[q] * W_ACI_F1_[fc][q][0][n]);
00562
00563 for (int q=0; q<nQuad(); q++, coeffPtr++){
00564 double f = (*coeffPtr)*detJ;
00565 for (int n=0; n<nNodes(); n++){
00566 aPtr[offset+n] += f*wi[n + nNodes()*q];
00567 }
00568 }
00569 }
00570 else
00571 {
00572 for (int q=0; q<nQuad(); q++, coeffPtr++){
00573 double f = (*coeffPtr)*detJ;
00574 for (int n=0; n<nNodes(); n++){
00575 aPtr[offset+n] += f*w[n + nNodes()*q];
00576 }
00577 }
00578 }
00579 }
00580 }
00581 else
00582 {
00583 for (int c=0; c<JVol.numCells(); c++, offset+=nNodes())
00584 {
00585 Tabs tab2;
00586 double detJ = fabs(JVol.detJ()[c]);
00587 int fc = 0;
00588 if (nFacetCases() != 1) fc = facetIndex[c];
00589 const Array<double>& w = W_[fc];
00590 SUNDANCE_MSG4(integrationVerb(), tab2 << "c=" << c << " detJ=" << detJ);
00591
00592 for (int q=0; q<nQuad(); q++, coeffPtr++)
00593 {
00594 Tabs tab3;
00595 double f = (*coeffPtr)*detJ;
00596 SUNDANCE_MSG4(integrationVerb(), tab3 << "q=" << q << " coeff=" <<
00597 *coeffPtr << " coeff*detJ=" << f);
00598 for (int n=0; n<nNodes(); n++)
00599 {
00600 Tabs tab4;
00601 SUNDANCE_MSG4(integrationVerb(), tab4 << "n=" << n << " w=" <<
00602 w[n + nNodes()*q]);
00603 aPtr[offset+n] += f*w[n + nNodes()*q];
00604 }
00605 }
00606 }
00607 }
00608 if (integrationVerb() >= 4)
00609 {
00610 Tabs tab2;
00611 Out::os() << tab2 << "integration results on cell:" << std::endl;
00612 Out::os() << tab2 << setw(10) << "n" << setw(12) << "I_n" << std::endl;
00613 for (int n=0; n<nNodes(); n++)
00614 {
00615 Out::os() << tab2 << setw(10) << n
00616 << setw(12) << setprecision(6) << aPtr[offset+n] << std::endl;
00617 }
00618 }
00619
00620 SUNDANCE_MSG5(integrationVerb(), tabs << "output A = ");
00621 if (integrationVerb() >= 5) writeTable(Out::os(), tabs, *A, 6);
00622 addFlops( JVol.numCells() * (1 + nQuad() * (1 + 2*nNodes())) );
00623 }
00624 else
00625 {
00626
00627
00628
00629
00630 createOneFormTransformationMatrix(JTrans, JVol);
00631
00632 SUNDANCE_MSG4(transformVerb(),
00633 Tabs() << "transformation matrix=" << G(alpha()));
00634
00635 double* GPtr = &(G(alpha())[0]);
00636
00637 if (useSumFirstMethod())
00638 {
00639 transformSummingFirst(JVol.numCells(), facetIndex, cellLIDs, GPtr, coeff, A);
00640 }
00641 else
00642 {
00643 transformSummingLast(JVol.numCells(), facetIndex, cellLIDs, GPtr, coeff, A);
00644 }
00645 }
00646 addFlops(flops);
00647 }
00648
00649
00650 void QuadratureIntegral::transformTwoForm(const CellJacobianBatch& JTrans,
00651 const CellJacobianBatch& JVol,
00652 const Array<int>& facetIndex,
00653 const RCP<Array<int> >& cellLIDs,
00654 const double* const coeff,
00655 RCP<Array<double> >& A) const
00656 {
00657 TimeMonitor timer(quadrature2Timer());
00658 Tabs tabs;
00659 TEUCHOS_TEST_FOR_EXCEPTION(order() != 2, std::logic_error,
00660 "QuadratureIntegral::transformTwoForm() called for form "
00661 "of order " << order());
00662 SUNDANCE_MSG2(integrationVerb(), tabs << "doing two form by quadrature");
00663
00664 const Array<int>* cellLID = cellLIDs.get();
00665
00666
00667
00668
00669 if (testDerivOrder() == 0 && unkDerivOrder() == 0)
00670 {
00671 double* aPtr = &((*A)[0]);
00672 double* coeffPtr = (double*) coeff;
00673 int offset = 0 ;
00674
00675 if (globalCurve().isCurveValid())
00676 {
00677
00678 Array<double> quadWeightsTmp = quadWeights_[0];
00679 Array<Point> quadPointsTmp = quadPts_[0];
00680 bool isCutByCurve;
00681
00682 for (int c=0; c<JVol.numCells(); c++, offset+=nNodes())
00683 {
00684 double detJ = fabs(JVol.detJ()[c]);
00685 int fc = 0;
00686 if (nFacetCases() != 1) fc = facetIndex[c];
00687 const Array<double>& w = W_[fc];
00688
00689 quadWeightsTmp = quadWeights_[fc];
00690 quadPointsTmp = quadPts_[fc];
00691
00692 quad_.getAdaptedWeights(cellType(), dim(), (*cellLID)[c] , fc ,
00693 mesh() , globalCurve() , quadPointsTmp , quadWeightsTmp , isCutByCurve );
00694 if (isCutByCurve){
00695 Array<double> wi;
00696 wi.resize(nQuad() * nNodesTest() *nNodesUnk() );
00697 for (int ii = 0 ; ii < wi.size() ; ii++ ) wi[ii] = 0.0;
00698 for (int nt = 0 ; nt < nNodesTest() ; nt++){
00699 for (int nu=0; nu<nNodesUnk(); nu++)
00700 for (int q=0 ; q < quadWeightsTmp.size() ; q++)
00701
00702 wi[nu + nNodesUnk()*(nt + nNodesTest()*q)] +=
00703 chop(quadWeightsTmp[q] * W_ACI_F2_[fc][q][0][nt][0][nu]);
00704 }
00705 for (int q=0; q<nQuad(); q++, coeffPtr++){
00706 double f = (*coeffPtr)*detJ;
00707 for (int n=0; n<nNodes(); n++){
00708 aPtr[offset+n] += f*wi[n + nNodes()*q];
00709 }
00710 }
00711 }
00712 else
00713 {
00714 for (int q=0; q<nQuad(); q++, coeffPtr++){
00715 double f = (*coeffPtr)*detJ;
00716 for (int n=0; n<nNodes(); n++){
00717 aPtr[offset+n] += f*w[n + nNodes()*q];
00718 }
00719 }
00720 }
00721 }
00722 }
00723 else
00724 {
00725
00726 for (int c=0; c<JVol.numCells(); c++, offset+=nNodes())
00727 {
00728 double detJ = fabs(JVol.detJ()[c]);
00729 int fc = 0;
00730 if (nFacetCases() != 1) fc = facetIndex[c];
00731 const Array<double>& w = W_[fc];
00732 for (int q=0; q<nQuad(); q++, coeffPtr++)
00733 {
00734 double f = (*coeffPtr)*detJ;
00735 for (int n=0; n<nNodes(); n++)
00736 {
00737 aPtr[offset+n] += f*w[n + nNodes()*q];
00738 }
00739 }
00740 }
00741 }
00742 addFlops( JVol.numCells() * (1 + nQuad() * (1 + 2*nNodes())) );
00743 }
00744 else
00745 {
00746 createTwoFormTransformationMatrix(JTrans, JVol);
00747 double* GPtr;
00748
00749 if (testDerivOrder() == 0)
00750 {
00751 GPtr = &(G(beta())[0]);
00752 SUNDANCE_MSG2(transformVerb(),
00753 Tabs() << "transformation matrix=" << G(beta()));
00754 }
00755 else if (unkDerivOrder() == 0)
00756 {
00757 GPtr = &(G(alpha())[0]);
00758 SUNDANCE_MSG2(transformVerb(),
00759 Tabs() << "transformation matrix=" << G(alpha()));
00760 }
00761 else
00762 {
00763 GPtr = &(G(alpha(), beta())[0]);
00764 SUNDANCE_MSG2(transformVerb(),
00765 Tabs() << "transformation matrix="
00766 << G(alpha(),beta()));
00767 }
00768
00769
00770 if (useSumFirstMethod())
00771 {
00772 transformSummingFirst(JTrans.numCells(), facetIndex, cellLIDs, GPtr, coeff, A);
00773 }
00774 else
00775 {
00776 transformSummingLast(JTrans.numCells(), facetIndex, cellLIDs, GPtr, coeff, A);
00777 }
00778 }
00779 }
00780
00781 void QuadratureIntegral
00782 ::transformSummingFirst(int nCells,
00783 const Array<int>& facetIndex,
00784 const RCP<Array<int> >& cellLIDs,
00785 const double* const GPtr,
00786 const double* const coeff,
00787 RCP<Array<double> >& A) const
00788 {
00789 double* aPtr = &((*A)[0]);
00790 double* coeffPtr = (double*) coeff;
00791 const Array<int>* cellLID = cellLIDs.get();
00792
00793 int transSize = 0;
00794 if (order()==2)
00795 {
00796 transSize = nRefDerivTest() * nRefDerivUnk();
00797 }
00798 else
00799 {
00800 transSize = nRefDerivTest();
00801 }
00802
00803
00804 static Array<double> sumWorkspace;
00805
00806 int swSize = transSize * nNodes();
00807 sumWorkspace.resize(swSize);
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819 if (globalCurve().isCurveValid())
00820 {
00821
00822 Array<double> quadWeightsTmp = quadWeights_[0];
00823 Array<Point> quadPointsTmp = quadPts_[0];
00824 bool isCutByCurve;
00825
00826 for (int c=0; c<nCells; c++)
00827 {
00828
00829 for (int i=0; i<swSize; i++) sumWorkspace[i]=0.0;
00830 int fc = 0;
00831 if (nFacetCases() > 1) fc = facetIndex[c];
00832
00833 const Array<double>& w = W_[fc];
00834
00835 quadWeightsTmp = quadWeights_[fc];
00836 quadPointsTmp = quadPts_[fc];
00837
00838 quad_.getAdaptedWeights(cellType(), dim(), (*cellLID)[c] , fc ,
00839 mesh() , globalCurve() , quadPointsTmp , quadWeightsTmp , isCutByCurve );
00840 if (isCutByCurve){
00841 Array<double> wi;
00842 if (order()==1){
00843 wi.resize(nQuad() * nRefDerivTest() * nNodesTest());
00844 for (int ii = 0 ; ii < wi.size() ; ii++ ) wi[ii] = 0.0;
00845 for (int n = 0 ; n < nNodes() ; n++)
00846 for (int t=0; t<nRefDerivTest(); t++)
00847 for (int q=0 ; q < quadWeightsTmp.size() ; q++)
00848
00849 wi[n + nNodes()*(t + nRefDerivTest()*q)] +=
00850 chop(quadWeightsTmp[q] * W_ACI_F1_[fc][q][t][n]);
00851 }else{
00852 wi.resize(nQuad() * nRefDerivTest() * nNodesTest()
00853 * nRefDerivUnk() * nNodesUnk() );
00854 for (int ii = 0 ; ii < wi.size() ; ii++ ) wi[ii] = 0.0;
00855
00856 for (int t=0; t<nRefDerivTest(); t++)
00857 for (int nt = 0 ; nt < nNodesTest() ; nt++){
00858 for (int u=0; u<nRefDerivUnk(); u++)
00859 for (int nu=0; nu<nNodesUnk(); nu++)
00860 for (int q=0 ; q < quadWeightsTmp.size() ; q++)
00861
00862
00863 wi[nu + nNodesUnk()*(nt + nNodesTest()*(u + nRefDerivUnk()*(t + nRefDerivTest()*q)))] +=
00864 chop(quadWeightsTmp[q] * W_ACI_F2_[fc][q][t][nt][u][nu]);
00865 }
00866 }
00867 for (int q=0; q<nQuad(); q++, coeffPtr++){
00868 double f = (*coeffPtr);
00869 for (int n=0; n<swSize; n++){
00870 sumWorkspace[n] += f*wi[n + q*swSize];
00871 }
00872 }
00873 }
00874 else{
00875 for (int q=0; q<nQuad(); q++, coeffPtr++)
00876 {
00877 double f = (*coeffPtr);
00878 for (int n=0; n<swSize; n++)
00879 {
00880 sumWorkspace[n] += f*w[n + q*swSize];
00881 }
00882 }
00883 }
00884 }
00885 }
00886 else
00887 {
00888 for (int c=0; c<nCells; c++)
00889 {
00890
00891 for (int i=0; i<swSize; i++) sumWorkspace[i]=0.0;
00892 int fc = 0;
00893 if (nFacetCases() > 1) fc = facetIndex[c];
00894
00895 const Array<double>& w = W_[fc];
00896
00897 for (int q=0; q<nQuad(); q++, coeffPtr++)
00898 {
00899 double f = (*coeffPtr);
00900 for (int n=0; n<swSize; n++)
00901 {
00902 sumWorkspace[n] += f*w[n + q*swSize];
00903 }
00904 }
00905
00906
00907 const double * const gCell = &(GPtr[transSize*c]);
00908 double* aCell = aPtr + nNodes()*c;
00909 for (int i=0; i<nNodes(); i++)
00910 {
00911 for (int j=0; j<transSize; j++)
00912 {
00913 aCell[i] += sumWorkspace[nNodes()*j + i] * gCell[j];
00914 }
00915 }
00916 }
00917 }
00918 int flops = 2*(nCells * nNodes() * transSize) * (nQuad() + 1) ;
00919 addFlops(flops);
00920 }
00921
00922 void QuadratureIntegral
00923 ::transformSummingLast(int nCells,
00924 const Array<int>& facetIndex,
00925 const RCP<Array<int> >& cellLIDs,
00926 const double* const GPtr,
00927 const double* const coeff,
00928 RCP<Array<double> >& A) const
00929 {
00930 double* aPtr = &((*A)[0]);
00931 int transSize = 0;
00932 const Array<int>* cellLID = cellLIDs.get();
00933
00934 if (order()==2)
00935 {
00936 transSize = nRefDerivTest() * nRefDerivUnk();
00937 }
00938 else
00939 {
00940 transSize = nRefDerivTest();
00941 }
00942
00943
00944
00945 static Array<double> jWorkspace;
00946 jWorkspace.resize(transSize);
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957 if (globalCurve().isCurveValid()){
00958
00959 Array<double> quadWeightsTmp = quadWeights_[0];
00960 Array<Point> quadPointsTmp = quadPts_[0];
00961 bool isCutByCurve;
00962
00963 for (int c=0; c<nCells; c++)
00964 {
00965 const double* const gCell = &(GPtr[transSize*c]);
00966 double* aCell = aPtr + nNodes()*c;
00967 int fc = 0;
00968 if (nFacetCases() > 1) fc = facetIndex[c];
00969 const Array<double>& w = W_[fc];
00970
00971 quadWeightsTmp = quadWeights_[fc];
00972 quadPointsTmp = quadPts_[fc];
00973
00974 quad_.getAdaptedWeights(cellType(), dim(), (*cellLID)[c] , fc ,
00975 mesh() , globalCurve() , quadPointsTmp , quadWeightsTmp , isCutByCurve );
00976 if (isCutByCurve){
00977 Array<double> wi;
00978 if (order()==1){
00979 wi.resize(nQuad() * nRefDerivTest() * nNodesTest());
00980 for (int ii = 0 ; ii < wi.size() ; ii++ ) wi[ii] = 0.0;
00981 for (int n = 0 ; n < nNodes() ; n++)
00982 for (int t=0; t<nRefDerivTest(); t++)
00983 for (int q=0 ; q < quadWeightsTmp.size() ; q++)
00984
00985 wi[n + nNodes()*(t + nRefDerivTest()*q)] +=
00986 chop(quadWeightsTmp[q] * W_ACI_F1_[fc][q][t][n]);
00987 }else{
00988 wi.resize(nQuad() * nRefDerivTest() * nNodesTest()
00989 * nRefDerivUnk() * nNodesUnk() );
00990 for (int ii = 0 ; ii < wi.size() ; ii++ ) wi[ii] = 0.0;
00991
00992 for (int t=0; t<nRefDerivTest(); t++)
00993 for (int nt = 0 ; nt < nNodesTest() ; nt++){
00994 for (int u=0; u<nRefDerivUnk(); u++)
00995 for (int nu=0; nu<nNodesUnk(); nu++)
00996 for (int q=0 ; q < quadWeightsTmp.size() ; q++)
00997
00998
00999 wi[nu + nNodesUnk()*(nt + nNodesTest()*(u + nRefDerivUnk()*(t + nRefDerivTest()*q)))] +=
01000 chop(quadWeightsTmp[q] * W_ACI_F2_[fc][q][t][nt][u][nu]);
01001 }
01002 }
01003 for (int q=0; q<nQuad(); q++){
01004 double f = coeff[c*nQuad() + q];
01005 for (int t=0; t<transSize; t++) jWorkspace[t]=f*gCell[t];
01006
01007 for (int n=0; n<nNodes(); n++){
01008 for (int t=0; t<transSize; t++){
01009 aCell[n] += jWorkspace[t]*wi[n + nNodes()*(t + transSize*q)];
01010 }
01011 }
01012 }
01013 }
01014 else{
01015 for (int q=0; q<nQuad(); q++){
01016 double f = coeff[c*nQuad() + q];
01017 for (int t=0; t<transSize; t++) jWorkspace[t]=f*gCell[t];
01018
01019 for (int n=0; n<nNodes(); n++){
01020 for (int t=0; t<transSize; t++){
01021 aCell[n] += jWorkspace[t]*w[n + nNodes()*(t + transSize*q)];
01022 }
01023 }
01024 }
01025 }
01026 }
01027 }
01028 else
01029 {
01030 for (int c=0; c<nCells; c++)
01031 {
01032 const double* const gCell = &(GPtr[transSize*c]);
01033 double* aCell = aPtr + nNodes()*c;
01034 int fc = 0;
01035 if (nFacetCases() > 1) fc = facetIndex[c];
01036 const Array<double>& w = W_[fc];
01037 for (int q=0; q<nQuad(); q++)
01038 {
01039 double f = coeff[c*nQuad() + q];
01040 for (int t=0; t<transSize; t++) jWorkspace[t]=f*gCell[t];
01041
01042 for (int n=0; n<nNodes(); n++)
01043 {
01044 for (int t=0; t<transSize; t++)
01045 {
01046 aCell[n] += jWorkspace[t]*w[n + nNodes()*(t + transSize*q)];
01047 }
01048 }
01049 }
01050 }
01051 }
01052 int flops = nCells * nQuad() * transSize * (1 + 2*nNodes()) ;
01053 addFlops(flops);
01054 }
