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 "SundanceQuadratureEvalMediator.hpp"
00032 #include "SundanceCoordExpr.hpp"
00033 #include "SundanceTempStack.hpp"
00034 #include "SundanceCellDiameterExpr.hpp"
00035 #include "SundanceCellVectorExpr.hpp"
00036 #include "SundanceSpatialDerivSpecifier.hpp"
00037 #include "SundanceDiscreteFunction.hpp"
00038 #include "SundanceDiscreteFuncElement.hpp"
00039 #include "SundanceCellJacobianBatch.hpp"
00040 #include "SundanceOut.hpp"
00041 #include "PlayaTabs.hpp"
00042 #include "PlayaExceptions.hpp"
00043
00044 #include "Teuchos_BLAS.hpp"
00045
00046
00047 using namespace Teuchos;
00048 TEUCHOS_TIMER(coordEvalTimer, "Quad mediator: coord eval")
00049
00050 using namespace Sundance;
00051 using namespace Playa;
00052 using std::endl;
00053 using std::setw;
00054
00055
00056 QuadratureEvalMediator
00057 ::QuadratureEvalMediator(const Mesh& mesh,
00058 int cellDim,
00059 const QuadratureFamily& quad)
00060 : StdFwkEvalMediator(mesh, cellDim),
00061 numEvaluationCases_(-1),
00062 quad_(quad),
00063 numQuadPtsForCellType_(),
00064 quadPtsForReferenceCell_(),
00065 quadPtsReferredToMaxCell_(),
00066 physQuadPts_(),
00067 refFacetBasisVals_(2)
00068 {}
00069
00070 Time& QuadratureEvalMediator::coordEvaluationTimer()
00071 {
00072 return coordEvalTimer();
00073 }
00074
00075 void QuadratureEvalMediator::setCellType(const CellType& cellType,
00076 const CellType& maxCellType,
00077 bool isInternBdry)
00078 {
00079 StdFwkEvalMediator::setCellType(cellType, maxCellType, isInternBdry);
00080
00081 Tabs tab;
00082
00083 SUNDANCE_MSG2(verb(), tab << "QuadEvalMed::setCellType: cellType="
00084 << cellType << " cellDim=" << cellDim() << " maxCellType=" << maxCellType);
00085 SUNDANCE_MSG2(verb(), tab << "integration spec: ="
00086 << integrationCellSpec());
00087 SUNDANCE_MSG2(verb(), tab << "forbid cofacet integrations ="
00088 << forbidCofacetIntegrations());
00089 if (isInternalBdry())
00090 {
00091 SUNDANCE_MSG2(verb(), tab << "working on internal boundary");
00092 }
00093
00094
00095
00096
00097 if (cellType != maxCellType)
00098 {
00099 Tabs tab1;
00100 SUNDANCE_MSG2(verb(), tab1 << "working out #facet cases");
00101 numEvaluationCases_ = numFacets(maxCellType, cellDim());
00102 }
00103 else
00104 {
00105 Tabs tab1;
00106 SUNDANCE_MSG2(verb(), tab1 << "no need for facet cases; work on original cell");
00107 numEvaluationCases_ = 1;
00108 }
00109 SUNDANCE_MSG2(verb(), tab << "num eval cases =" << numEvaluationCases_);
00110
00111 if (!isInternalBdry()
00112 && quadPtsReferredToMaxCell_.containsKey(cellType)) return;
00113
00114 if (!quadPtsForReferenceCell_.containsKey(cellType))
00115 {
00116 SUNDANCE_MSG2(verb(), tab << "creating quad points for ref cell type="
00117 << cellType);
00118 RCP<Array<Point> > pts = rcp(new Array<Point>());
00119 RCP<Array<double> > wgts = rcp(new Array<double>());
00120
00121 quad_.getPoints(cellType, *pts, *wgts);
00122 quadPtsForReferenceCell_.put(cellType, pts);
00123
00124 numQuadPtsForCellType_.put(cellType, pts->size());
00125 }
00126
00127 if (!quadPtsReferredToMaxCell_.containsKey(cellType))
00128 {
00129 SUNDANCE_MSG2(verb(), tab << "creating quad points for max cell type="
00130 << maxCellType);
00131 RCP<Array<Array<Point> > > facetPts
00132 = rcp(new Array<Array<Point> >(numEvaluationCases()));
00133 RCP<Array<Array<double> > > facetWgts
00134 = rcp(new Array<Array<double> >(numEvaluationCases()));
00135
00136 for (int fc=0; fc<numEvaluationCases(); fc++)
00137 {
00138 if (cellType != maxCellType)
00139 {
00140 quad_.getFacetPoints(maxCellType, cellDim(), fc,
00141 (*facetPts)[fc], (*facetWgts)[fc]);
00142 }
00143 else
00144 {
00145 quad_.getPoints(maxCellType, (*facetPts)[fc], (*facetWgts)[fc]);
00146 }
00147 }
00148 quadPtsReferredToMaxCell_.put(cellType, facetPts);
00149 }
00150 }
00151
00152 int QuadratureEvalMediator::numQuadPts(const CellType& ct) const
00153 {
00154 TEUCHOS_TEST_FOR_EXCEPTION(!numQuadPtsForCellType_.containsKey(ct),
00155 std::runtime_error,
00156 "no quadrature points have been tabulated for cell type=" << ct);
00157 return numQuadPtsForCellType_.get(ct);
00158 }
00159
00160 void QuadratureEvalMediator::evalCellDiameterExpr(const CellDiameterExpr* expr,
00161 RCP<EvalVector>& vec) const
00162 {
00163 Tabs tabs;
00164 SUNDANCE_MSG2(verb(),tabs
00165 << "QuadratureEvalMediator evaluating cell diameter expr "
00166 << expr->toString());
00167
00168 int nQuad = numQuadPts(cellType());
00169 int nCells = cellLID()->size();
00170
00171 SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad);
00172 Array<double> diameters;
00173 mesh().getCellDiameters(cellDim(), *cellLID(), diameters);
00174
00175 vec->resize(nQuad*nCells);
00176 double * const xx = vec->start();
00177 int k=0;
00178 for (int c=0; c<nCells; c++)
00179 {
00180 double h = diameters[c];
00181 for (int q=0; q<nQuad; q++, k++)
00182 {
00183 xx[k] = h;
00184 }
00185 }
00186 }
00187
00188 void QuadratureEvalMediator::evalCellVectorExpr(const CellVectorExpr* expr,
00189 RCP<EvalVector>& vec) const
00190 {
00191 Tabs tabs;
00192 SUNDANCE_MSG2(verb(),tabs
00193 << "QuadratureEvalMediator evaluating cell vector expr "
00194 << expr->toString());
00195
00196 int nQuad = numQuadPts(cellType());
00197 int nCells = cellLID()->size();
00198
00199 vec->resize(nQuad*nCells);
00200
00201 SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad);
00202 int dir = expr->componentIndex();
00203
00204 Array<Point> vectors;
00205 if (expr->isNormal())
00206 {
00207 mesh().outwardNormals(*cellLID(), vectors);
00208 }
00209 else
00210 {
00211 TEUCHOS_TEST_FOR_EXCEPTION(cellDim() != 1, std::runtime_error,
00212 "unable to compute tangent vectors for cell dim = " << cellDim());
00213 mesh().tangentsToEdges(*cellLID(), vectors);
00214 }
00215
00216 double * const xx = vec->start();
00217 int k=0;
00218 for (int c=0; c<nCells; c++)
00219 {
00220 double n = vectors[c][dir];
00221 for (int q=0; q<nQuad; q++, k++)
00222 {
00223 xx[k] = n;
00224 }
00225 }
00226 }
00227
00228 void QuadratureEvalMediator::evalCoordExpr(const CoordExpr* expr,
00229 RCP<EvalVector>& vec) const
00230 {
00231 Tabs tabs;
00232 SUNDANCE_MSG2(verb(),tabs
00233 << "QuadratureEvalMediator evaluating coord expr "
00234 << expr->toString());
00235
00236 TimeMonitor timer(coordEvalTimer());
00237
00238 computePhysQuadPts();
00239 int nQuad = physQuadPts_.length();
00240 int d = expr->dir();
00241
00242 SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad);
00243
00244 vec->resize(nQuad);
00245 double * const xx = vec->start();
00246 for (int q=0; q<nQuad; q++)
00247 {
00248 xx[q] = physQuadPts_[q][d];
00249 }
00250 }
00251
00252 RCP<Array<Array<Array<double> > > > QuadratureEvalMediator
00253 ::getFacetRefBasisVals(const BasisFamily& basis, int diffOrder) const
00254 {
00255 Tabs tab;
00256 RCP<Array<Array<Array<double> > > > rtn ;
00257
00258 CellType evalCellType = cellType();
00259 if (cellDim() != maxCellDim())
00260 {
00261 if (verb() >= 2)
00262 {
00263 Out::os() << tab << "alwaysUseCofacets = "
00264 << ElementIntegral::alwaysUseCofacets() << std::endl;
00265 Out::os() << tab << "diffOrder = " << diffOrder << std::endl;
00266 }
00267 if (ElementIntegral::alwaysUseCofacets() || diffOrder>0)
00268 {
00269 if (!cofacetCellsAreReady()) setupFacetTransformations();
00270 evalCellType = maxCellType();
00271
00272 TEUCHOS_TEST_FOR_EXCEPTION(!cofacetCellsAreReady(), std::runtime_error,
00273 "cofacet cells not ready in getFacetRefBasisVals()");
00274 }
00275 }
00276
00277 SUNDANCE_MSG2(verb(), tab << "eval cell type = " << evalCellType);
00278 typedef OrderedPair<BasisFamily, CellType> key;
00279
00280 int nDerivResults = 1;
00281 if (diffOrder==1) nDerivResults = maxCellDim();
00282
00283 if (!refFacetBasisVals_[diffOrder].containsKey(key(basis, cellType())))
00284 {
00285 SUNDANCE_OUT(this->verb() > 2,
00286 tab << "computing basis values on facet quad pts");
00287 rtn = rcp(new Array<Array<Array<double> > >(numEvaluationCases()));
00288
00289 Array<Array<Array<Array<double> > > > tmp(nDerivResults);
00290
00291 if (verb() >= 2)
00292 {
00293 Out::os() << tab << "numEvalCases = " << numEvaluationCases()
00294 << std::endl;
00295 Out::os() << tab << "diff order = " << diffOrder << std::endl;
00296 Out::os() << tab << "cell type = " << cellType() << std::endl;
00297 Out::os() << tab << "quad pt map = ";
00298 if (evalCellType!=cellType())
00299 {
00300 Out::os() << quadPtsReferredToMaxCell_ << std::endl;
00301 }
00302 else
00303 {
00304 Out::os() << quadPtsForReferenceCell_ << std::endl;
00305 }
00306 }
00307
00308 if (evalCellType!=cellType())
00309 {
00310 TEUCHOS_TEST_FOR_EXCEPTION(quadPtsReferredToMaxCell_.size() == 0,
00311 std::runtime_error,
00312 "empty quadrature point map (max cell)");
00313 }
00314 else
00315 {
00316 TEUCHOS_TEST_FOR_EXCEPTION(quadPtsForReferenceCell_.size() == 0,
00317 std::runtime_error,
00318 "empty quadrature point map (ref cell)");
00319 }
00320
00321 for (int fc=0; fc<numEvaluationCases(); fc++)
00322 {
00323 Tabs tab1;
00324 (*rtn)[fc].resize(basis.dim());
00325 SUNDANCE_MSG2(verb(), tab1 << "fc = " << fc);
00326
00327 for (int r=0; r<nDerivResults; r++)
00328 {
00329 tmp[r].resize(basis.dim());
00330 MultiIndex mi;
00331 if (diffOrder==1)
00332 {
00333 mi[r]=1;
00334 }
00335 SpatialDerivSpecifier deriv(mi);
00336
00337
00338 if (evalCellType != cellType())
00339 {
00340 SUNDANCE_MSG2(verb(), tab1 << "referring to max cell");
00341 basis.refEval(evalCellType,
00342 (*(quadPtsReferredToMaxCell_.get(cellType())))[fc],
00343 deriv, tmp[r], verb());
00344 }
00345 else
00346 {
00347 SUNDANCE_MSG2(verb(), tab1 << "computing on reference cell");
00348 basis.refEval(evalCellType,
00349 (*(quadPtsForReferenceCell_.get(cellType()))),
00350 deriv, tmp[r], verb());
00351 }
00352 }
00353
00354
00355
00356 int dim = maxCellDim();
00357 int nQuad = tmp[0][0].size();
00358 int nNodes = tmp[0][0][0].size();
00359 int nTot = dim * nQuad * nNodes;
00360 for (int d=0; d<basis.dim(); d++)
00361 {
00362 (*rtn)[fc][d].resize(nTot);
00363 for (int r=0; r<nDerivResults; r++)
00364 {
00365 for (int q=0; q<nQuad; q++)
00366 {
00367 for (int n=0; n<nNodes; n++)
00368 {
00369 (*rtn)[fc][d][(n*nQuad + q)*nDerivResults + r] = tmp[r][d][q][n];
00370 }
00371 }
00372 }
00373 }
00374 }
00375 refFacetBasisVals_[diffOrder].put(key(basis, cellType()), rtn);
00376 }
00377 else
00378 {
00379 SUNDANCE_OUT(this->verb() > 2,
00380 tab << "reusing facet basis values on quad pts");
00381 rtn = refFacetBasisVals_[diffOrder].get(key(basis, cellType()));
00382 }
00383
00384 return rtn;
00385 }
00386
00387 Array<Array<double> >* QuadratureEvalMediator
00388 ::getRefBasisVals(const BasisFamily& basis, int diffOrder) const
00389 {
00390 Tabs tab;
00391 RCP<Array<Array<Array<double> > > > fRtn
00392 = getFacetRefBasisVals(basis, diffOrder);
00393 return &((*fRtn)[0]);
00394 }
00395
00396
00397 void QuadratureEvalMediator
00398 ::evalDiscreteFuncElement(const DiscreteFuncElement* expr,
00399 const Array<MultiIndex>& multiIndices,
00400 Array<RCP<EvalVector> >& vec) const
00401 {
00402 const DiscreteFunctionData* f = DiscreteFunctionData::getData(expr);
00403 TEUCHOS_TEST_FOR_EXCEPTION(f==0, std::logic_error,
00404 "QuadratureEvalMediator::evalDiscreteFuncElement() called "
00405 "with expr that is not a discrete function");
00406 Tabs tab;
00407
00408 SUNDANCE_MSG1(verb(),tab << "QuadEvalMed evaluating DF " << expr->name());
00409
00410 int nQuad = numQuadPts(cellType());
00411 int myIndex = expr->myIndex();
00412
00413 for (int i=0; i<multiIndices.size(); i++)
00414 {
00415 Tabs tab1;
00416 const MultiIndex& mi = multiIndices[i];
00417 SUNDANCE_MSG2(dfVerb(),
00418 tab1 << "evaluating DF " << expr->name()
00419 << " for multiindex " << mi << std::endl
00420 << tab1 << "num cells = " << cellLID()->size() << std::endl
00421 << tab1 << "num quad points = " << nQuad << std::endl
00422 << tab1 << "my index = " << expr->myIndex() << std::endl
00423 << tab1 << "num funcs = " << f->discreteSpace().nFunc());
00424
00425 vec[i]->resize(cellLID()->size() * nQuad);
00426
00427 if (mi.order() == 0)
00428 {
00429 Tabs tab2;
00430 SUNDANCE_MSG2(dfVerb(),tab2 << "in mi.order() == 0 branch");
00431 if (!fCache().containsKey(f) || !fCacheIsValid()[f])
00432 {
00433 fillFunctionCache(f, mi);
00434 }
00435 else
00436 {
00437 SUNDANCE_MSG2(dfVerb(),tab2 << "reusing function cache");
00438 }
00439
00440 const RCP<const MapStructure>& mapStruct = mapStructCache()[f];
00441 int chunk = mapStruct->chunkForFuncID(myIndex);
00442 int funcIndex = mapStruct->indexForFuncID(myIndex);
00443 int nFuncs = mapStruct->numFuncs(chunk);
00444
00445 SUNDANCE_MSG3(dfVerb(),tab2 << "chunk number = " << chunk << std::endl
00446 << tab2 << "function index=" << funcIndex << " of nFuncs="
00447 << nFuncs);
00448
00449 const RCP<Array<Array<double> > >& cacheVals
00450 = fCache()[f];
00451
00452 SUNDANCE_MSG4(dfVerb(),tab2 << "cached function values=" << (*cacheVals)[chunk]);
00453
00454 const double* cachePtr = &((*cacheVals)[chunk][0]);
00455 double* vecPtr = vec[i]->start();
00456
00457 int cellSize = nQuad*nFuncs;
00458 int offset = funcIndex*nQuad;
00459 SUNDANCE_MSG3(dfVerb(),tab2 << "cell size=" << cellSize << ", offset="
00460 << offset);
00461 int k = 0;
00462 for (int c=0; c<cellLID()->size(); c++)
00463 {
00464 for (int q=0; q<nQuad; q++, k++)
00465 {
00466 vecPtr[k] = cachePtr[c*cellSize + offset + q];
00467 }
00468 }
00469 SUNDANCE_MSG4(dfVerb(),tab2 << "result vector=");
00470 if (dfVerb() >= 5)
00471 {
00472 vec[i]->print(Out::os());
00473 computePhysQuadPts();
00474 k=0;
00475 for (int c=0; c<cellLID()->size(); c++)
00476 {
00477 Out::os() << tab2 << "-------------------------------------------"
00478 << std::endl;
00479 Out::os() << tab2 << "c=" << c << " cell LID=" << (*cellLID())[c]
00480 << std::endl;
00481 Tabs tab3;
00482 for (int q=0; q<nQuad; q++, k++)
00483 {
00484 Out::os() << tab3 << "q=" << q << " " << physQuadPts_[k]
00485 << " val=" << vecPtr[k] << std::endl;
00486 }
00487 }
00488 }
00489 }
00490 else
00491 {
00492 Tabs tab2;
00493 SUNDANCE_MSG2(dfVerb(),tab2 << "in mi.order() != 0 branch");
00494 if (!dfCache().containsKey(f) || !dfCacheIsValid()[f])
00495 {
00496 fillFunctionCache(f, mi);
00497 }
00498 else
00499 {
00500 SUNDANCE_MSG3(dfVerb(),tab2 << "reusing function cache");
00501 }
00502
00503 RCP<const MapStructure> mapStruct = mapStructCache()[f];
00504 int chunk = mapStruct->chunkForFuncID(myIndex);
00505 int funcIndex = mapStruct->indexForFuncID(myIndex);
00506 int nFuncs = mapStruct->numFuncs(chunk);
00507
00508
00509 SUNDANCE_MSG3(dfVerb(),tab2 << "chunk number = " << chunk << std::endl
00510 << tab2 << "function index=" << funcIndex << " of nFuncs="
00511 << nFuncs);
00512
00513 const RCP<Array<Array<double> > >& cacheVals
00514 = dfCache()[f];
00515
00516 SUNDANCE_MSG4(dfVerb(),tab2 << "cached function values=" << (*cacheVals)[chunk]);
00517
00518 int dim = maxCellDim();
00519 int pDir = mi.firstOrderDirection();
00520 const double* cachePtr = &((*cacheVals)[chunk][0]);
00521 double* vecPtr = vec[i]->start();
00522
00523 int cellSize = nQuad*nFuncs*dim;
00524 int offset = funcIndex * nQuad * dim;
00525 int k = 0;
00526
00527 SUNDANCE_MSG2(dfVerb(),tab2 << "dim=" << dim << ", pDir=" << pDir
00528 << ", cell size=" << cellSize << ", offset="
00529 << offset);
00530 for (int c=0; c<cellLID()->size(); c++)
00531 {
00532 for (int q=0; q<nQuad; q++, k++)
00533 {
00534 vecPtr[k] = cachePtr[c*cellSize + offset + q*dim + pDir];
00535 }
00536 }
00537 SUNDANCE_MSG4(dfVerb(),tab2 << "result vector=");
00538 if (dfVerb() >= 5)
00539 {
00540 vec[i]->print(Out::os());
00541 computePhysQuadPts();
00542 k=0;
00543 for (int c=0; c<cellLID()->size(); c++)
00544 {
00545 Out::os() << tab2 << "-------------------------------------------"
00546 << std::endl;
00547 Out::os() << tab2 << "c=" << c << " cell LID=" << (*cellLID())[c]
00548 << std::endl;
00549 Tabs tab3;
00550 for (int q=0; q<nQuad; q++, k++)
00551 {
00552 Out::os() << tab3 << "q=" << q << " " << physQuadPts_[k]
00553 << " val=" << vecPtr[k] << std::endl;
00554 }
00555 }
00556 }
00557 }
00558 }
00559 }
00560
00561 void QuadratureEvalMediator::fillFunctionCache(const DiscreteFunctionData* f,
00562 const MultiIndex& mi) const
00563 {
00564 Tabs tab0(0);
00565
00566 SUNDANCE_MSG2(dfVerb(), tab0 << "QuadratureEvalMediator::fillFunctionCache()");
00567 SUNDANCE_MSG2(dfVerb(), tab0 << "multiIndex=" << mi);
00568
00569
00570 int diffOrder = mi.order();
00571 CellType evalCellType = cellType();
00572
00573 int funcSupportDim = f->map()->cellDim();
00574
00575 TEUCHOS_TEST_FOR_EXCEPT(diffOrder > 0 && funcSupportDim < maxCellDim());
00576
00577 int flops = 0;
00578 double jFlops = CellJacobianBatch::totalFlops();
00579
00580 RCP<Array<Array<double> > > localValues;
00581 RCP<const MapStructure> mapStruct;
00582
00583 Teuchos::BLAS<int,double> blas;
00584
00585 {
00586 Tabs tab1;
00587 if (cellDim() != maxCellDim())
00588 {
00589 if (dfVerb() >= 2)
00590 {
00591 Out::os() << tab1 << "alwaysUseCofacets = "
00592 << ElementIntegral::alwaysUseCofacets() << std::endl;
00593 Out::os() << tab1 << "diffOrder = " << diffOrder << std::endl;
00594 }
00595 if (diffOrder==0 && funcSupportDim < maxCellDim())
00596 {
00597 evalCellType = cellType();
00598 }
00599 else
00600 {
00601 evalCellType = maxCellType();
00602 }
00603 }
00604
00605 SUNDANCE_MSG2(dfVerb(), tab1 << "cell type=" << cellType());
00606 SUNDANCE_MSG2(dfVerb(), tab1 << "max cell type=" << maxCellType());
00607 SUNDANCE_MSG2(dfVerb(), tab1 << "eval cell type=" << evalCellType);
00608
00609 SUNDANCE_MSG2(dfVerb(), tab1 << "packing local values");
00610
00611 if (!localValueCacheIsValid().containsKey(f)
00612 || !localValueCacheIsValid().get(f))
00613 {
00614 Tabs tab2;
00615 SUNDANCE_MSG2(dfVerb(), tab2 << "filling cache");
00616 localValues = rcp(new Array<Array<double> >());
00617 if (cellDim() != maxCellDim())
00618 {
00619 if (diffOrder==0 && funcSupportDim < maxCellDim())
00620 {
00621 mapStruct = f->getLocalValues(cellDim(), *cellLID(),
00622 *localValues);
00623 }
00624 else
00625 {
00626 if (!cofacetCellsAreReady()) setupFacetTransformations();
00627 mapStruct = f->getLocalValues(maxCellDim(), *cofacetCellLID(),
00628 *localValues);
00629 }
00630 }
00631 else
00632 {
00633 mapStruct = f->getLocalValues(maxCellDim(), *cellLID(), *localValues);
00634 }
00635
00636 TEUCHOS_TEST_FOR_EXCEPT(mapStruct.get() == 0);
00637 localValueCache().put(f, localValues);
00638 mapStructCache().put(f, mapStruct);
00639 localValueCacheIsValid().put(f, true);
00640 }
00641 else
00642 {
00643 Tabs tab2;
00644 SUNDANCE_MSG2(dfVerb(), tab2 << "reusing local value cache");
00645 localValues = localValueCache().get(f);
00646 mapStruct = mapStructCache().get(f);
00647 }
00648 }
00649
00650 RCP<Array<Array<double> > > cacheVals;
00651
00652 if (mi.order()==0)
00653 {
00654 if (fCache().containsKey(f))
00655 {
00656 cacheVals = fCache().get(f);
00657 }
00658 else
00659 {
00660 cacheVals = rcp(new Array<Array<double> >(mapStruct->numBasisChunks()));
00661 fCache().put(f, cacheVals);
00662 }
00663 fCacheIsValid().put(f, true);
00664 }
00665 else
00666 {
00667 if (dfCache().containsKey(f))
00668 {
00669 cacheVals = dfCache().get(f);
00670 }
00671 else
00672 {
00673 cacheVals = rcp(new Array<Array<double> >(mapStruct->numBasisChunks()));
00674 dfCache().put(f, cacheVals);
00675 }
00676 dfCacheIsValid().put(f, true);
00677 }
00678
00679
00680
00681 for (int chunk=0; chunk<mapStruct->numBasisChunks(); chunk++)
00682 {
00683 Tabs tab1;
00684 SUNDANCE_MSG2(dfVerb(), tab1 << "processing dof map chunk=" << chunk
00685 << " of " << mapStruct->numBasisChunks());
00686 Tabs tab2;
00687 BasisFamily basis = rcp_dynamic_cast<BasisFamilyBase>(mapStruct->basis(chunk));
00688 SUNDANCE_MSG4(dfVerb(), tab2 << "basis=" << basis);
00689
00690 int nFuncs = mapStruct->numFuncs(chunk);
00691 SUNDANCE_MSG2(dfVerb(), tab2 << "num funcs in this chunk=" << nFuncs);
00692
00693
00694 Array<double>& cache = (*cacheVals)[chunk];
00695
00696 int nQuad = numQuadPts(cellType());
00697 int nCells = cellLID()->size();
00698 SUNDANCE_MSG2(dfVerb(), tab2 << "num quad points=" << nQuad);
00699 SUNDANCE_MSG2(dfVerb(), tab2 << "num cells=" << nCells);
00700
00701
00702 int nDir;
00703
00704 if (mi.order()==1)
00705 {
00706 nDir = maxCellDim();
00707 }
00708 else
00709 {
00710 nDir = 1;
00711 }
00712 cache.resize(cellLID()->size() * nQuad * nDir * nFuncs);
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741 if (cellType() != evalCellType)
00742 {
00743 Tabs tab2;
00744 SUNDANCE_MSG2(dfVerb(),
00745 tab2 << "evaluating by reference to maximal cell");
00746
00747 RCP<Array<Array<Array<double> > > > refFacetBasisValues
00748 = getFacetRefBasisVals(basis, mi.order());
00749
00750
00751
00752
00753
00754
00755 int nNodes = basis.nReferenceDOFsWithFacets(maxCellType(), maxCellType());
00756 int nRowsA = nQuad*nDir;
00757 int nColsA = nNodes;
00758 int nColsB = nFuncs;
00759 int lda = nRowsA;
00760 int ldb = nNodes;
00761 int ldc = lda;
00762 double alpha = 1.0;
00763 double beta = 0.0;
00764
00765 SUNDANCE_MSG2(dfVerb(), tab2 << "building transformation matrices cell-by-cell");
00766 int vecComp = 0;
00767 for (int c=0; c<nCells; c++)
00768 {
00769 int facetIndex = facetIndices()[c];
00770 double* A = &((*refFacetBasisValues)[facetIndex][vecComp][0]);
00771 double* B = &((*localValues)[chunk][c*nNodes*nFuncs]);
00772 double* C = &((*cacheVals)[chunk][c*nRowsA*nColsB]);
00773 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, nRowsA, nColsB, nColsA,
00774 alpha, A, lda, B, ldb, beta, C, ldc);
00775 }
00776 }
00777 else
00778 {
00779
00780
00781
00782
00783 Tabs tab2;
00784 SUNDANCE_MSG2(dfVerb(), tab2 << "building batch transformation matrix");
00785
00786 Array<Array<double> >* refBasisValues
00787 = getRefBasisVals(basis, diffOrder);
00788 int nNodes = basis.nReferenceDOFsWithFacets(maxCellType(), cellType());
00789 int nRowsA = nQuad*nDir;
00790 int nColsA = nNodes;
00791 int nColsB = nFuncs*nCells;
00792 int lda = nRowsA;
00793 int ldb = nNodes;
00794 int ldc = lda;
00795 double alpha = 1.0;
00796 double beta = 0.0;
00797 int vecComp = 0;
00798 if (dfVerb() >= 5)
00799 {
00800 Tabs tab3;
00801 Out::os() << tab2 << "Printing values at nodes" << std::endl;
00802 for (int c=0; c<nCells; c++)
00803 {
00804 Out::os() << tab3 << "-------------------------------------------"
00805 << std::endl;
00806 Out::os() << tab3 << "c=" << c << " cell LID=" << (*cellLID())[c]
00807 << std::endl;
00808 Tabs tab4;
00809 int offset = c*nNodes*nFuncs;
00810 for (int n=0; n<nNodes; n++)
00811 {
00812 Out::os() << tab4 << "n=" << n;
00813 for (int f=0; f<nFuncs; f++)
00814 {
00815 Out::os() << " " << (*localValues)[chunk][offset + n*nFuncs + f];
00816 }
00817 Out::os() << std::endl;
00818 }
00819 }
00820 Out::os() << tab2 << "-------------------------------------------";
00821 Out::os() << tab2 << "Printing reference basis at nodes" << std::endl;
00822 Out::os() << tab2 << "-------------------------------------------"
00823 << std::endl;
00824 for (int n=0; n<nNodes; n++)
00825 {
00826 Out::os() << tab3 << "node=" << n << std::endl;
00827 for (int q=0; q<nQuad; q++)
00828 {
00829 Tabs tab4;
00830 Out::os() << tab4 << "q=" << q;
00831 for (int r=0; r<nDir; r++)
00832 {
00833 Out::os () << " "
00834 << (*refBasisValues)[vecComp][(n*nQuad + q)*nDir + r];
00835 }
00836 Out::os() << std::endl;
00837 }
00838 }
00839 }
00840 double* A = &((*refBasisValues)[vecComp][0]);
00841 double* B = &((*localValues)[chunk][0]);
00842 double* C = &((*cacheVals)[chunk][0]);
00843 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, nRowsA, nColsB, nColsA, alpha,
00844 A, lda, B, ldb, beta, C, ldc );
00845 }
00846
00847
00848
00849 const CellJacobianBatch& J = JTrans();
00850 double* C = &((*cacheVals)[chunk][0]);
00851 if (mi.order()==1)
00852 {
00853 SUNDANCE_MSG2(dfVerb(), tab1 << "doing transformations via dgemm");
00854 Tabs tab2;
00855 SUNDANCE_MSG2(dfVerb(), tab2 << "Jacobian batch nCells=" << J.numCells());
00856 SUNDANCE_MSG2(dfVerb(), tab2 << "Jacobian batch cell dim=" << J.cellDim());
00857 SUNDANCE_MSG2(dfVerb(), tab2 << "Jacobian batch spatial dim=" << J.spatialDim());
00858
00859 int nRhs = nQuad * nFuncs;
00860 for (int c=0; c<cellLID()->size(); c++)
00861 {
00862 double* rhsPtr = &(C[(nRhs * nDir)*c]);
00863 J.applyInvJ(c, 0, rhsPtr, nRhs, false);
00864 }
00865 }
00866 else
00867 {
00868 SUNDANCE_MSG2(dfVerb(), tab1 << "no derivatives to transform");
00869 }
00870 SUNDANCE_MSG2(dfVerb(), tab1 << "done transformations");
00871 }
00872
00873 jFlops = CellJacobianBatch::totalFlops() - jFlops;
00874 addFlops(flops + jFlops);
00875
00876 SUNDANCE_MSG2(dfVerb(),
00877 tab0 << "done QuadratureEvalMediator::fillFunctionCache()");
00878 }
00879
00880 void QuadratureEvalMediator::computePhysQuadPts() const
00881 {
00882 if (cacheIsValid())
00883 {
00884 Tabs tab0(0);
00885 SUNDANCE_MSG2(verb(), tab0 << "reusing phys quad points");
00886 }
00887 else
00888 {
00889 Tabs tab0(0);
00890 double jFlops = CellJacobianBatch::totalFlops();
00891 SUNDANCE_MSG2(verb(), tab0 << "computing phys quad points");
00892 physQuadPts_.resize(0);
00893 if (cellDim() != maxCellDim() && ElementIntegral::alwaysUseCofacets()
00894 && !isInternalBdry())
00895 {
00896 Tabs tab1;
00897 SUNDANCE_MSG2(verb(), tab1 << "using cofacets");
00898 SUNDANCE_MSG3(verb(), tab1 << "num cofacets = " << cofacetCellLID()->size());
00899 Array<Point> tmp;
00900 Array<int> cell(1);
00901 for (int c=0; c<cellLID()->size(); c++)
00902 {
00903 cell[0] = (*cofacetCellLID())[c];
00904 int facetIndex = facetIndices()[c];
00905 const Array<Point>& refFacetPts
00906 = (*(quadPtsReferredToMaxCell_.get(cellType())))[facetIndex];
00907 tmp.resize(refFacetPts.size());
00908 mesh().pushForward(maxCellDim(), cell, refFacetPts, tmp);
00909 for (int q=0; q<tmp.size(); q++)
00910 {
00911 physQuadPts_.append(tmp[q]);
00912 }
00913 }
00914 }
00915 else
00916 {
00917 const Array<Point>& refPts
00918 = *(quadPtsForReferenceCell_.get(cellType()));
00919 mesh().pushForward(cellDim(), *cellLID(),
00920 refPts, physQuadPts_);
00921 }
00922 addFlops(CellJacobianBatch::totalFlops() - jFlops);
00923 cacheIsValid() = true;
00924 SUNDANCE_OUT(this->verb() > 2,
00925 "phys quad: " << physQuadPts_);
00926 }
00927 }
00928
00929
00930 void QuadratureEvalMediator::print(std::ostream& os) const
00931 {
00932 if (cacheIsValid())
00933 {
00934 Tabs tab0;
00935 os << tab0 << "Physical quadrature points" << std::endl;
00936 int i=0;
00937 for (int c=0; c<cellLID()->size(); c++)
00938 {
00939 Tabs tab1;
00940 os << tab1 << "cell " << c << std::endl;
00941 for (int q=0; q<physQuadPts_.size()/cellLID()->size(); q++, i++)
00942 {
00943 Tabs tab2;
00944 os << tab2 << "q=" << q << " " << physQuadPts_[i] << std::endl;
00945 }
00946 }
00947 }
00948 }
00949
00950
00951 void QuadratureEvalMediator
00952 ::showResults(std::ostream& os,
00953 const RCP<SparsitySuperset>& sp,
00954 const Array<RCP<EvalVector> >& vecResults,
00955 const Array<double>& constantResults) const
00956 {
00957 Tabs tabs(0);
00958
00959
00960
00961
00962
00963 int maxlen = 25;
00964 for (int i=0; i<sp->numDerivs(); i++)
00965 {
00966 int s = sp->deriv(i).toString().length();
00967 if (s > maxlen) maxlen = s;
00968 }
00969
00970
00971 int vecIndex=0;
00972 int constIndex = 0;
00973 os << tabs << "Results Superset" << std::endl;
00974 for (int i=0; i<sp->numDerivs(); i++)
00975 {
00976 Tabs tab1;
00977 os << tab1 << i << " " ;
00978 os.width(maxlen);
00979 os.setf(std::ios_base::left, std::ios_base::adjustfield);
00980 os << sp->deriv(i).toString() ;
00981 Tabs tab2;
00982 switch(sp->state(i))
00983 {
00984 case ZeroDeriv:
00985 os << "Zero" << std::endl;
00986 break;
00987 case ConstantDeriv:
00988 os << "const val=" << constantResults[constIndex++] << std::endl;
00989 break;
00990 case VectorDeriv:
00991 if (vecResults[vecIndex].get()==0)
00992 {
00993 os << "{Null}";
00994 }
00995 else
00996 {
00997 const double* data = vecResults[vecIndex]->start();
00998 if (EvalVector::shadowOps())
00999 {
01000 TEUCHOS_TEST_FOR_EXCEPT(vecResults[vecIndex]->str().size()==0);
01001 os << vecResults[vecIndex]->str();
01002 }
01003 os << endl;
01004 int nQuad = numQuadPts(cellType());
01005 int k=0;
01006 TEUCHOS_TEST_FOR_EXCEPT(cellLID()->size() * nQuad
01007 != vecResults[vecIndex]->length());
01008 for (int c=0; c<cellLID()->size(); c++)
01009 {
01010 Tabs tab3;
01011 os << tab3 << "cell LID=" << (*cellLID())[c] << endl;
01012 for (int q=0; q<nQuad; q++, k++)
01013 {
01014 Tabs tab4;
01015 os << tab4 << "q=" << setw(5) << q
01016 << setw(10) << Utils::chop(data[k]) << endl;
01017 }
01018 }
01019 }
01020 vecIndex++;
01021 os << std::endl;
01022 break;
01023 }
01024 }
01025 }