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