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