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