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 "SundanceReducedIntegral.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& reduced0IntegrationTimer()
00059 {
00060 static RCP<Time> rtn
00061 = TimeMonitor::getNewTimer("reduced 0-form integration");
00062 return *rtn;
00063 }
00064
00065
00066 static Time& reduced1IntegrationTimer()
00067 {
00068 static RCP<Time> rtn
00069 = TimeMonitor::getNewTimer("reduced 1-form integration");
00070 return *rtn;
00071 }
00072
00073
00074 static Time& reduced2IntegrationTimer()
00075 {
00076 static RCP<Time> rtn
00077 = TimeMonitor::getNewTimer("reduced 2-form integration");
00078 return *rtn;
00079 }
00080
00081
00082
00083 ReducedIntegral::ReducedIntegral(int spatialDim,
00084 const CellType& maxCellType,
00085 int dim,
00086 const CellType& cellType,
00087 const QuadratureFamily& quad_in,
00088 bool isInternalBdry,
00089 const ParametrizedCurve& globalCurve,
00090 const Mesh& mesh,
00091 int verb)
00092 : ElementIntegral(spatialDim, maxCellType, dim, cellType, isInternalBdry, globalCurve , mesh,
00093 verb), W_()
00094 {
00095 Tabs tab0(0);
00096
00097 SUNDANCE_MSG1(setupVerb(),
00098 tab0 << "************* creating reduced 0-form integrals ********");
00099 if (setupVerb()) describe(Out::os());
00100
00101
00102
00103 QuadratureFamily quad = new GaussianQuadrature(1);
00104
00105
00106 TEUCHOS_TEST_FOR_EXCEPT(globalCurve.isCurveValid());
00107
00108 Array<Point> quadPts;
00109 Array<double> quadWeights;
00110
00111 W_.resize(1);
00112 W_[0].resize(1);
00113
00114 quad.getPoints(cellType, quadPts, quadWeights);
00115
00116 for (int q=0; q<quadWeights.size(); q++) {
00117 W_[0][0] += quadWeights[q];
00118 }
00119 }
00120
00121 ReducedIntegral::ReducedIntegral(int spatialDim,
00122 const CellType& maxCellType,
00123 int dim,
00124 const CellType& cellType,
00125 const BasisFamily& testBasis,
00126 int alpha,
00127 int testDerivOrder,
00128 const QuadratureFamily& quad_in,
00129 bool isInternalBdry,
00130 const ParametrizedCurve& globalCurve,
00131 const Mesh& mesh,
00132 int verb)
00133 : ElementIntegral(spatialDim, maxCellType, dim, cellType,
00134 testBasis, alpha, testDerivOrder, isInternalBdry, globalCurve , mesh , verb), W_()
00135 {
00136 Tabs tab0(0);
00137 SUNDANCE_MSG1(setupVerb(),
00138 tab0 << "************* creating reduced 1-form integrals ********");
00139 if (setupVerb()) describe(Out::os());
00140 assertLinearForm();
00141
00142 W_.resize(nFacetCases());
00143
00144
00145 QuadratureType qType = new GaussianQuadratureType();
00146 int reqOrder = qType.findValidOrder(cellType,
00147 std::max(1, testBasis.order()));
00148
00149 SUNDANCE_MSG2(setupVerb(),
00150 tab0 << "using quadrature order=" << reqOrder);
00151
00152
00153 QuadratureFamily quad = qType.createQuadFamily(reqOrder);
00154
00155
00156 TEUCHOS_TEST_FOR_EXCEPT(globalCurve.isCurveValid());
00157
00158
00159
00160
00161 for (int fc=0; fc<nFacetCases(); fc++)
00162 {
00163 Tabs tab1;
00164 SUNDANCE_MSG2(setupVerb(),
00165 tab1 << "evaluation case=" << fc << " of " << nFacetCases());
00166
00167 W_[fc].resize(nRefDerivTest() * nNodesTest());
00168
00169
00170 for (int i=0; i<W_[fc].size(); i++) { W_[fc][i]=0.0; }
00171
00172 Array<Array<Array<Array<double> > > > testBasisVals(nRefDerivTest());
00173
00174
00175
00176 Array<Point> quadPts;
00177 Array<double> quadWeights;
00178 getQuad(quad, fc, quadPts, quadWeights);
00179
00180 int nQuad = quadPts.size();
00181
00182
00183 for (int r=0; r<nRefDerivTest(); r++)
00184 {
00185 Tabs tab2;
00186 SUNDANCE_MSG2(setupVerb(), tab2 << "evaluating basis derivative "
00187 << r << " of " << nRefDerivTest());
00188
00189 testBasisVals[r].resize(testBasis.dim());
00190 MultiIndex mi;
00191 if (testDerivOrder==1) mi[r] = 1;
00192 SpatialDerivSpecifier deriv(mi);
00193 testBasis.refEval(evalCellType(), quadPts, deriv,
00194 testBasisVals[r], setupVerb());
00195 }
00196
00197
00198 SUNDANCE_MSG2(setupVerb(), tab1 << "doing quadrature");
00199 int vecComp = 0;
00200 for (int q=0; q<nQuad; q++)
00201 {
00202 for (int t=0; t<nRefDerivTest(); t++)
00203 {
00204 for (int nt=0; nt<nNodesTest(); nt++)
00205 {
00206 value(fc, t, nt)
00207 += chop(quadWeights[q] * testBasisVals[t][vecComp][q][nt]) ;
00208 }
00209 }
00210 }
00211
00212 for (int i=0; i<W_[fc].size(); i++) W_[fc][i] = chop(W_[fc][i]);
00213
00214 addFlops(3*nQuad*nRefDerivTest()*nNodesTest() + W_[fc].size());
00215 }
00216
00217
00218 SUNDANCE_MSG4(setupVerb(), tab0 << "--------------------------------------");
00219 SUNDANCE_MSG4(setupVerb(), tab0 << "reduced linear form integral results");
00220 if (setupVerb() >= 4)
00221 {
00222 for (int fc=0; fc<nFacetCases(); fc++)
00223 {
00224 Tabs tab1;
00225 SUNDANCE_MSG4(setupVerb(), tab1 << "------ evaluation case " << fc << " of "
00226 << nFacetCases() << "-------");
00227
00228 for (int r=0; r<nRefDerivTest(); r++)
00229 {
00230 Tabs tab2;
00231
00232 MultiIndex mi;
00233 if (testDerivOrder==1) mi[r] = 1;
00234 SUNDANCE_MSG1(setupVerb(), tab2 << "multiindex=" << mi);
00235
00236 ios_base::fmtflags oldFlags = Out::os().flags();
00237 Out::os().setf(ios_base::right);
00238 Out::os().setf(ios_base::showpoint);
00239 for (int nt=0; nt<nNodesTest(); nt++)
00240 {
00241 Tabs tab3;
00242 Out::os() << tab3 << setw(10) << nt
00243 << setw(12) << std::setprecision(5) << value(fc, r, nt)
00244 << std::endl;
00245 }
00246 Out::os().flags(oldFlags);
00247 }
00248 }
00249 }
00250
00251 SUNDANCE_MSG1(setupVerb(), tab0 << "done reduced linear form ctor");
00252 }
00253
00254
00255
00256
00257 ReducedIntegral::ReducedIntegral(int spatialDim,
00258 const CellType& maxCellType,
00259 int dim,
00260 const CellType& cellType,
00261 const BasisFamily& testBasis,
00262 int alpha,
00263 int testDerivOrder,
00264 const BasisFamily& unkBasis,
00265 int beta,
00266 int unkDerivOrder,
00267 const QuadratureFamily& quad_in,
00268 bool isInternalBdry,
00269 const ParametrizedCurve& globalCurve,
00270 const Mesh& mesh,
00271 int verb)
00272 : ElementIntegral(spatialDim, maxCellType, dim, cellType,
00273 testBasis, alpha, testDerivOrder,
00274 unkBasis, beta, unkDerivOrder, isInternalBdry, globalCurve , mesh ,verb), W_()
00275
00276 {
00277 Tabs tab0(0);
00278 SUNDANCE_MSG1(setupVerb(),
00279 tab0 << "************* creating reduced 2-form integrals ********");
00280 if (setupVerb()) describe(Out::os());
00281
00282 assertBilinearForm();
00283
00284 W_.resize(nFacetCases());
00285
00286 QuadratureType qType = new GaussianQuadratureType();
00287 int reqOrder = qType.findValidOrder(cellType,
00288 std::max(1, unkBasis.order() + testBasis.order()));
00289
00290 SUNDANCE_MSG2(setupVerb(),
00291 tab0 << "using quadrature order=" << reqOrder);
00292 QuadratureFamily quad = qType.createQuadFamily(reqOrder);
00293
00294
00295 TEUCHOS_TEST_FOR_EXCEPT(globalCurve.isCurveValid());
00296
00297
00298 SUNDANCE_MSG2(setupVerb(),
00299 tab0 << "processing evaluation cases");
00300
00301 for (int fc=0; fc<nFacetCases(); fc++)
00302 {
00303 Tabs tab1;
00304 SUNDANCE_MSG1(setupVerb(), tab1 << "------ evaluation case " << fc << " of "
00305 << nFacetCases() << "-------");
00306
00307 W_[fc].resize(nRefDerivTest() * nNodesTest() * nRefDerivUnk() * nNodesUnk());
00308 for (int i=0; i<W_[fc].size(); i++) W_[fc][i]=0.0;
00309
00310 Array<Array<Array<Array<double> > > > testBasisVals(nRefDerivTest());
00311 Array<Array<Array<Array<double> > > > unkBasisVals(nRefDerivUnk());
00312
00313 Array<Point> quadPts;
00314 Array<double> quadWeights;
00315 getQuad(quad, fc, quadPts, quadWeights);
00316
00317 int nQuad = quadPts.size();
00318
00319 for (int r=0; r<nRefDerivTest(); r++)
00320 {
00321 Tabs tab2;
00322 SUNDANCE_MSG2(setupVerb(), tab2
00323 << "evaluating test function basis derivative "
00324 << r << " of " << nRefDerivTest());
00325 testBasisVals[r].resize(testBasis.dim());
00326 MultiIndex mi;
00327 if (testDerivOrder==1) mi[r] = 1;
00328 SpatialDerivSpecifier deriv(mi);
00329 testBasis.refEval(evalCellType(), quadPts, deriv,
00330 testBasisVals[r], setupVerb());
00331 }
00332
00333 for (int r=0; r<nRefDerivUnk(); r++)
00334 {
00335 Tabs tab2;
00336 SUNDANCE_MSG2(setupVerb(), tab2
00337 << "evaluating unknown function basis derivative "
00338 << r << " of " << nRefDerivUnk());
00339 unkBasisVals[r].resize(unkBasis.dim());
00340 MultiIndex mi;
00341 if (unkDerivOrder==1) mi[r] = 1;
00342 SpatialDerivSpecifier deriv(mi);
00343 unkBasis.refEval(evalCellType(),
00344 quadPts, deriv, unkBasisVals[r], setupVerb());
00345 }
00346
00347 SUNDANCE_MSG2(setupVerb(), tab1 << "doing quadrature...");
00348 int vecComp = 0;
00349 for (int q=0; q<nQuad; q++)
00350 {
00351 for (int t=0; t<nRefDerivTest(); t++)
00352 {
00353 for (int nt=0; nt<nNodesTest(); nt++)
00354 {
00355 for (int u=0; u<nRefDerivUnk(); u++)
00356 {
00357 for (int nu=0; nu<nNodesUnk(); nu++)
00358 {
00359 value(fc, t, nt, u, nu)
00360 += chop(quadWeights[q] * testBasisVals[t][vecComp][q][nt]
00361 * unkBasisVals[u][vecComp][q][nu]);
00362 }
00363 }
00364 }
00365 }
00366 }
00367 SUNDANCE_MSG2(setupVerb(), tab1 << "...done");
00368 addFlops(4*nQuad*nRefDerivTest()*nNodesTest()*nRefDerivUnk()*nNodesUnk()
00369 + W_[fc].size());
00370 for (int i=0; i<W_[fc].size(); i++) W_[fc][i] = chop(W_[fc][i]);
00371 }
00372
00373 SUNDANCE_MSG1(setupVerb(), tab0
00374 << "----------------------------------------");
00375 SUNDANCE_MSG4(setupVerb(), tab0
00376 << "reduced bilinear form integral results");
00377 if (setupVerb() >= 4)
00378 {
00379 for (int fc=0; fc<nFacetCases(); fc++)
00380 {
00381 Tabs tab1;
00382 SUNDANCE_MSG4(setupVerb(), tab1 << "evaluation case " << fc << " of "
00383 << nFacetCases());
00384
00385 for (int rt=0; rt<nRefDerivTest(); rt++)
00386 {
00387 for (int ru=0; ru<nRefDerivUnk(); ru++)
00388 {
00389 Tabs tab2;
00390 MultiIndex miTest;
00391 if (testDerivOrder==1) miTest[rt] = 1;
00392 MultiIndex miUnk;
00393 if (unkDerivOrder==1) miUnk[ru] = 1;
00394 SUNDANCE_MSG1(setupVerb(), tab2 << "test multiindex=" << miTest
00395 << " unk multiindex=" << miUnk);
00396
00397 ios_base::fmtflags oldFlags = Out::os().flags();
00398 Out::os().setf(ios_base::right);
00399 Out::os().setf(ios_base::showpoint);
00400 for (int nt=0; nt<nNodesTest(); nt++)
00401 {
00402 Tabs tab3;
00403 Out::os() << tab3 << setw(10) << nt;
00404 for (int nu=0; nu<nNodesUnk(); nu++)
00405 {
00406 Out::os() << setw(12) << std::setprecision(5)
00407 << value(fc, rt, nt, ru, nu) ;
00408 }
00409 Out::os() << std::endl;
00410 }
00411 Out::os().flags(oldFlags);
00412 }
00413 }
00414 }
00415 }
00416
00417 SUNDANCE_MSG1(setupVerb(), tab0 << "done reduced bilinear form ctor");
00418 }
00419
00420
00421
00422
00423 void ReducedIntegral::transformZeroForm(const CellJacobianBatch& JTrans,
00424 const CellJacobianBatch& JVol,
00425 const Array<int>& isLocalFlag,
00426 const Array<int>& facetIndex,
00427 const RCP<Array<int> >& cellLIDs,
00428 const double* const coeffs,
00429 RCP<Array<double> >& A) const
00430 {
00431 TimeMonitor timer(reduced0IntegrationTimer());
00432
00433 TEUCHOS_TEST_FOR_EXCEPTION(order() != 0, std::logic_error,
00434 "ReducedIntegral::transformZeroForm() called "
00435 "for form of order " << order());
00436
00437 Tabs tabs;
00438 SUNDANCE_MSG1(integrationVerb(), tabs << "doing zero form by reduced integration");
00439
00440 double& a = (*A)[0];
00441 int flops = 0;
00442
00443
00444
00445
00446 double w = W_[0][0];
00447 if ((int) isLocalFlag.size()==0)
00448 {
00449
00450 if (globalCurve().isCurveValid())
00451 {
00452 TEUCHOS_TEST_FOR_EXCEPT(globalCurve().isCurveValid());
00453 }
00454 else
00455 {
00456 for (int c=0; c<JVol.numCells(); c++)
00457 {
00458 flops+=3;
00459 a += w * coeffs[c] * fabs(JVol.detJ()[c]);
00460 }
00461 }
00462
00463 }
00464 else
00465 {
00466 TEUCHOS_TEST_FOR_EXCEPTION( (int) isLocalFlag.size() != JVol.numCells(),
00467 std::runtime_error,
00468 "mismatch between isLocalFlag.size()="
00469 << isLocalFlag.size()
00470 << " and JVol.numCells()=" << JVol.numCells());
00471
00472 if (globalCurve().isCurveValid())
00473 {
00474 TEUCHOS_TEST_FOR_EXCEPT(globalCurve().isCurveValid());
00475 }
00476 else
00477 {
00478 for (int c=0; c<JVol.numCells(); c++)
00479 {
00480 if (isLocalFlag[c])
00481 {
00482 flops+=3;
00483 a += w * coeffs[c] * fabs(JVol.detJ()[c]);
00484 }
00485 }
00486 }
00487 }
00488 addFlops(flops);
00489 }
00490
00491 void ReducedIntegral::transformOneForm(const CellJacobianBatch& JTrans,
00492 const CellJacobianBatch& JVol,
00493 const Array<int>& facetIndex,
00494 const RCP<Array<int> >& cellLIDs,
00495 const double* const coeffs,
00496 RCP<Array<double> >& A) const
00497 {
00498 TimeMonitor timer(reduced1IntegrationTimer());
00499 TEUCHOS_TEST_FOR_EXCEPTION(order() != 1, std::logic_error,
00500 "ReducedIntegral::transformOneForm() called for form "
00501 "of order " << order());
00502 Tabs tabs;
00503 SUNDANCE_MSG1(integrationVerb(),
00504 tabs << "doing one form by reduced integration");
00505
00506 int nfc = nFacetCases();
00507
00508
00509
00510
00511 if (testDerivOrder() == 0)
00512 {
00513 double* aPtr = &((*A)[0]);
00514 int count = 0;
00515 if (globalCurve().isCurveValid())
00516 {
00517 TEUCHOS_TEST_FOR_EXCEPT(globalCurve().isCurveValid());
00518 }
00519 else
00520 {
00521 for (int c=0; c<JVol.numCells(); c++)
00522 {
00523 double w_cell = coeffs[c] * fabs(JVol.detJ()[c]);
00524 int fc = 0;
00525 if (nfc != 1) fc = facetIndex[c];
00526
00527 const Array<double>& w = W_[fc];
00528 for (int n=0; n<nNodes(); n++, count++)
00529 {
00530 aPtr[count] += w_cell*w[n];
00531 }
00532 }
00533 }
00534 addFlops(JVol.numCells() * (nNodes() + 2));
00535 }
00536 else
00537 {
00538
00539
00540
00541 int nCells = JVol.numCells();
00542 double one = 1.0;
00543 int nTransRows = nRefDerivTest();
00544
00545 createOneFormTransformationMatrix(JTrans, JVol);
00546
00547 SUNDANCE_MSG3(transformVerb(),
00548 Tabs() << "transformation matrix=" << G(alpha()));
00549 int nNodes0 = nNodes();
00550
00551 if (nFacetCases()==1)
00552 {
00553
00554
00555
00556 if (globalCurve().isCurveValid())
00557 {
00558 TEUCHOS_TEST_FOR_EXCEPT(globalCurve().isCurveValid());
00559 }
00560 else
00561 {
00562 int N = 1;
00563 for (int c=0; c<JVol.numCells(); c++)
00564 {
00565 double* aPtr = &((*A)[c*nNodes0]);
00566 ::dgemm_("N", "N", &nNodes0, &N, &nTransRows, &(coeffs[c]), &(W_[0][0]),
00567 &nNodes0, &(G(alpha())[c*nTransRows]), &nTransRows, &one,
00568 aPtr, &nNodes0);
00569 }
00570 }
00571 }
00572 else
00573 {
00574
00575
00576 if (globalCurve().isCurveValid())
00577 {
00578 TEUCHOS_TEST_FOR_EXCEPT(globalCurve().isCurveValid());
00579 }
00580 else
00581 {
00582 int N = 1;
00583 for (int c=0; c<JVol.numCells(); c++)
00584 {
00585 int fc = 0;
00586 if (nfc != 1) fc = facetIndex[c];
00587 double* aPtr = &((*A)[c*nNodes0]);
00588 ::dgemm_("N", "N", &nNodes0, &N, &nTransRows, &(coeffs[c]), &(W_[fc][0]),
00589 &nNodes0, &(G(alpha())[c*nTransRows]), &nTransRows, &one,
00590 aPtr, &nNodes0);
00591 }
00592 }
00593 }
00594
00595 addFlops(2 * nNodes0 * nCells * nTransRows);
00596 }
00597 }
00598
00599 void ReducedIntegral::transformTwoForm(const CellJacobianBatch& JTrans,
00600 const CellJacobianBatch& JVol,
00601 const Array<int>& facetIndex,
00602 const RCP<Array<int> >& cellLIDs,
00603 const double* const coeffs,
00604 RCP<Array<double> >& A) const
00605 {
00606 TimeMonitor timer(reduced2IntegrationTimer());
00607 TEUCHOS_TEST_FOR_EXCEPTION(order() != 2, std::logic_error,
00608 "ReducedIntegral::transformTwoForm() called for form "
00609 "of order " << order());
00610
00611 Tabs tabs;
00612 SUNDANCE_MSG1(transformVerb(), tabs << "doing two form by reduced integration");
00613
00614 int nfc = nFacetCases();
00615
00616
00617
00618 if (testDerivOrder() == 0 && unkDerivOrder() == 0)
00619 {
00620 if (globalCurve().isCurveValid())
00621 {
00622 TEUCHOS_TEST_FOR_EXCEPT(globalCurve().isCurveValid());
00623 }
00624 else
00625 {
00626 double* aPtr = &((*A)[0]);
00627 int count = 0;
00628 for (int c=0; c<JVol.numCells(); c++)
00629 {
00630 int fc = 0;
00631 if (nFacetCases() != 1) fc = facetIndex[c];
00632
00633 const Array<double>& w = W_[fc];
00634 double w_cell = coeffs[c] * fabs(JVol.detJ()[c]);
00635 for (int n=0; n<nNodes(); n++, count++)
00636 {
00637 aPtr[count] += w_cell*w[n];
00638 }
00639 }
00640 }
00641 addFlops(JVol.numCells() * (nNodes() + 1));
00642 }
00643 else
00644 {
00645
00646
00647
00648 int nCells = JVol.numCells();
00649 double one = 1.0;
00650 int nTransRows = nRefDerivUnk()*nRefDerivTest();
00651
00652 createTwoFormTransformationMatrix(JTrans, JVol);
00653
00654 double* GPtr;
00655 if (testDerivOrder() == 0)
00656 {
00657 GPtr = &(G(beta())[0]);
00658 SUNDANCE_MSG2(transformVerb(),
00659 Tabs() << "transformation matrix=" << G(beta()));
00660 }
00661 else if (unkDerivOrder() == 0)
00662 {
00663 GPtr = &(G(alpha())[0]);
00664 SUNDANCE_MSG2(transformVerb(),
00665 Tabs() << "transformation matrix=" << G(alpha()));
00666 }
00667 else
00668 {
00669 GPtr = &(G(alpha(), beta())[0]);
00670 SUNDANCE_MSG2(transformVerb(),
00671 Tabs() << "transformation matrix="
00672 << G(alpha(),beta()));
00673 }
00674
00675 int nNodes0 = nNodes();
00676
00677 if (nFacetCases()==1)
00678 {
00679
00680
00681
00682 if (globalCurve().isCurveValid())
00683 {
00684 TEUCHOS_TEST_FOR_EXCEPT(globalCurve().isCurveValid());
00685 }
00686 else
00687 {
00688 int N = 1;
00689 for (int c=0; c<JVol.numCells(); c++)
00690 {
00691 double* aPtr = &((*A)[c*nNodes0]);
00692 double* gPtr = &(GPtr[c*nTransRows]);
00693 SUNDANCE_MSG2(integrationVerb(),
00694 tabs << "transforming c=" << c << ", W=" << W_[0]);
00695
00696 ::dgemm_("N", "N", &nNodes0, &N, &nTransRows, &(coeffs[c]), &(W_[0][0]),
00697 &nNodes0, gPtr, &nTransRows, &one,
00698 aPtr, &nNodes0);
00699 }
00700 }
00701 }
00702 else
00703 {
00704
00705
00706 if (globalCurve().isCurveValid())
00707 {
00708 TEUCHOS_TEST_FOR_EXCEPT(globalCurve().isCurveValid());
00709 }
00710 else
00711 {
00712 int N = 1;
00713 for (int c=0; c<JVol.numCells(); c++)
00714 {
00715 int fc = 0;
00716 if (nfc != 1) fc = facetIndex[c];
00717 double* aPtr = &((*A)[c*nNodes0]);
00718 double* gPtr = &(GPtr[c*nTransRows]);
00719 SUNDANCE_MSG2(integrationVerb(),
00720 tabs << "c=" << c << ", facet case=" << fc
00721 << " W=" << W_[fc]);
00722
00723 ::dgemm_("N", "N", &nNodes0, &N, &nTransRows, &(coeffs[c]), &(W_[fc][0]),
00724 &nNodes0, gPtr, &nTransRows, &one,
00725 aPtr, &nNodes0);
00726 }
00727 }
00728 }
00729
00730 addFlops(2 * nNodes0 * nCells * nTransRows);
00731 }
00732 }
00733