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 flops = 0;
00574 double jFlops = CellJacobianBatch::totalFlops();
00575
00576 RCP<Array<Array<double> > > localValues;
00577 RCP<const MapStructure> mapStruct;
00578
00579 Teuchos::BLAS<int,double> blas;
00580
00581 {
00582 Tabs tab1;
00583 if (cellDim() != maxCellDim())
00584 {
00585 if (dfVerb() >= 2)
00586 {
00587 Out::os() << tab1 << "alwaysUseCofacets = "
00588 << ElementIntegral::alwaysUseCofacets() << std::endl;
00589 Out::os() << tab1 << "diffOrder = " << diffOrder << std::endl;
00590 }
00591 if (ElementIntegral::alwaysUseCofacets() || diffOrder>0)
00592 {
00593 evalCellType = maxCellType();
00594 }
00595 }
00596
00597 SUNDANCE_MSG2(dfVerb(), tab1 << "cell type=" << cellType());
00598 SUNDANCE_MSG2(dfVerb(), tab1 << "max cell type=" << maxCellType());
00599 SUNDANCE_MSG2(dfVerb(), tab1 << "eval cell type=" << evalCellType);
00600
00601 SUNDANCE_MSG2(dfVerb(), tab1 << "packing local values");
00602
00603 if (!localValueCacheIsValid().containsKey(f)
00604 || !localValueCacheIsValid().get(f))
00605 {
00606 Tabs tab2;
00607 SUNDANCE_MSG2(dfVerb(), tab2 << "filling cache");
00608 localValues = rcp(new Array<Array<double> >());
00609 if (cellDim() != maxCellDim())
00610 {
00611 if (ElementIntegral::alwaysUseCofacets() || diffOrder>0)
00612 {
00613 if (!cofacetCellsAreReady()) setupFacetTransformations();
00614 mapStruct = f->getLocalValues(maxCellDim(), *cofacetCellLID(),
00615 *localValues);
00616 }
00617 else
00618 {
00619 mapStruct = f->getLocalValues(cellDim(), *cellLID(),
00620 *localValues);
00621 }
00622 }
00623 else
00624 {
00625 mapStruct = f->getLocalValues(maxCellDim(), *cellLID(), *localValues);
00626 }
00627
00628 TEUCHOS_TEST_FOR_EXCEPT(mapStruct.get() == 0);
00629 localValueCache().put(f, localValues);
00630 mapStructCache().put(f, mapStruct);
00631 localValueCacheIsValid().put(f, true);
00632 }
00633 else
00634 {
00635 Tabs tab2;
00636 SUNDANCE_MSG2(dfVerb(), tab2 << "reusing local value cache");
00637 localValues = localValueCache().get(f);
00638 mapStruct = mapStructCache().get(f);
00639 }
00640 }
00641
00642 RCP<Array<Array<double> > > cacheVals;
00643
00644 if (mi.order()==0)
00645 {
00646 if (fCache().containsKey(f))
00647 {
00648 cacheVals = fCache().get(f);
00649 }
00650 else
00651 {
00652 cacheVals = rcp(new Array<Array<double> >(mapStruct->numBasisChunks()));
00653 fCache().put(f, cacheVals);
00654 }
00655 fCacheIsValid().put(f, true);
00656 }
00657 else
00658 {
00659 if (dfCache().containsKey(f))
00660 {
00661 cacheVals = dfCache().get(f);
00662 }
00663 else
00664 {
00665 cacheVals = rcp(new Array<Array<double> >(mapStruct->numBasisChunks()));
00666 dfCache().put(f, cacheVals);
00667 }
00668 dfCacheIsValid().put(f, true);
00669 }
00670
00671
00672
00673 for (int chunk=0; chunk<mapStruct->numBasisChunks(); chunk++)
00674 {
00675 Tabs tab1;
00676 SUNDANCE_MSG2(dfVerb(), tab1 << "processing dof map chunk=" << chunk
00677 << " of " << mapStruct->numBasisChunks());
00678 Tabs tab2;
00679 BasisFamily basis = rcp_dynamic_cast<BasisFamilyBase>(mapStruct->basis(chunk));
00680 SUNDANCE_MSG4(dfVerb(), tab2 << "basis=" << basis);
00681
00682 int nFuncs = mapStruct->numFuncs(chunk);
00683 SUNDANCE_MSG2(dfVerb(), tab2 << "num funcs in this chunk=" << nFuncs);
00684
00685
00686 Array<double>& cache = (*cacheVals)[chunk];
00687
00688 int nQuad = numQuadPts(cellType());
00689 int nCells = cellLID()->size();
00690 SUNDANCE_MSG2(dfVerb(), tab2 << "num quad points=" << nQuad);
00691 SUNDANCE_MSG2(dfVerb(), tab2 << "num cells=" << nCells);
00692
00693
00694 int nDir;
00695
00696 if (mi.order()==1)
00697 {
00698 nDir = maxCellDim();
00699 }
00700 else
00701 {
00702 nDir = 1;
00703 }
00704 cache.resize(cellLID()->size() * nQuad * nDir * nFuncs);
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733 if (cellType() != evalCellType)
00734 {
00735 Tabs tab2;
00736 SUNDANCE_MSG2(dfVerb(),
00737 tab2 << "evaluating by reference to maximal cell");
00738
00739 RCP<Array<Array<Array<double> > > > refFacetBasisValues
00740 = getFacetRefBasisVals(basis, mi.order());
00741
00742
00743
00744
00745
00746
00747 int nNodes = basis.nReferenceDOFsWithFacets(maxCellType(), maxCellType());
00748 int nRowsA = nQuad*nDir;
00749 int nColsA = nNodes;
00750 int nColsB = nFuncs;
00751 int lda = nRowsA;
00752 int ldb = nNodes;
00753 int ldc = lda;
00754 double alpha = 1.0;
00755 double beta = 0.0;
00756
00757 SUNDANCE_MSG2(dfVerb(), tab2 << "building transformation matrices cell-by-cell");
00758 int vecComp = 0;
00759 for (int c=0; c<nCells; c++)
00760 {
00761 int facetIndex = facetIndices()[c];
00762 double* A = &((*refFacetBasisValues)[facetIndex][vecComp][0]);
00763 double* B = &((*localValues)[chunk][c*nNodes*nFuncs]);
00764 double* C = &((*cacheVals)[chunk][c*nRowsA*nColsB]);
00765 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, nRowsA, nColsB, nColsA,
00766 alpha, A, lda, B, ldb, beta, C, ldc);
00767 }
00768 }
00769 else
00770 {
00771
00772
00773
00774
00775 Tabs tab2;
00776 SUNDANCE_MSG2(dfVerb(), tab2 << "building batch transformation matrix");
00777
00778 Array<Array<double> >* refBasisValues
00779 = getRefBasisVals(basis, diffOrder);
00780 int nNodes = basis.nReferenceDOFsWithFacets(maxCellType(), cellType());
00781 int nRowsA = nQuad*nDir;
00782 int nColsA = nNodes;
00783 int nColsB = nFuncs*nCells;
00784 int lda = nRowsA;
00785 int ldb = nNodes;
00786 int ldc = lda;
00787 double alpha = 1.0;
00788 double beta = 0.0;
00789 int vecComp = 0;
00790 if (dfVerb() >= 5)
00791 {
00792 Tabs tab3;
00793 Out::os() << tab2 << "Printing values at nodes" << std::endl;
00794 for (int c=0; c<nCells; c++)
00795 {
00796 Out::os() << tab3 << "-------------------------------------------"
00797 << std::endl;
00798 Out::os() << tab3 << "c=" << c << " cell LID=" << (*cellLID())[c]
00799 << std::endl;
00800 Tabs tab4;
00801 int offset = c*nNodes*nFuncs;
00802 for (int n=0; n<nNodes; n++)
00803 {
00804 Out::os() << tab4 << "n=" << n;
00805 for (int f=0; f<nFuncs; f++)
00806 {
00807 Out::os() << " " << (*localValues)[chunk][offset + n*nFuncs + f];
00808 }
00809 Out::os() << std::endl;
00810 }
00811 }
00812 Out::os() << tab2 << "-------------------------------------------";
00813 Out::os() << tab2 << "Printing reference basis at nodes" << std::endl;
00814 Out::os() << tab2 << "-------------------------------------------"
00815 << std::endl;
00816 for (int n=0; n<nNodes; n++)
00817 {
00818 Out::os() << tab3 << "node=" << n << std::endl;
00819 for (int q=0; q<nQuad; q++)
00820 {
00821 Tabs tab4;
00822 Out::os() << tab4 << "q=" << q;
00823 for (int r=0; r<nDir; r++)
00824 {
00825 Out::os () << " "
00826 << (*refBasisValues)[vecComp][(n*nQuad + q)*nDir + r];
00827 }
00828 Out::os() << std::endl;
00829 }
00830 }
00831 }
00832 double* A = &((*refBasisValues)[vecComp][0]);
00833 double* B = &((*localValues)[chunk][0]);
00834 double* C = &((*cacheVals)[chunk][0]);
00835 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, nRowsA, nColsB, nColsA, alpha,
00836 A, lda, B, ldb, beta, C, ldc );
00837 }
00838
00839
00840
00841 const CellJacobianBatch& J = JTrans();
00842 double* C = &((*cacheVals)[chunk][0]);
00843 if (mi.order()==1)
00844 {
00845 SUNDANCE_MSG2(dfVerb(), tab1 << "doing transformations via dgemm");
00846 Tabs tab2;
00847 SUNDANCE_MSG2(dfVerb(), tab2 << "Jacobian batch nCells=" << J.numCells());
00848 SUNDANCE_MSG2(dfVerb(), tab2 << "Jacobian batch cell dim=" << J.cellDim());
00849 SUNDANCE_MSG2(dfVerb(), tab2 << "Jacobian batch spatial dim=" << J.spatialDim());
00850
00851 int nRhs = nQuad * nFuncs;
00852 for (int c=0; c<cellLID()->size(); c++)
00853 {
00854 double* rhsPtr = &(C[(nRhs * nDir)*c]);
00855 J.applyInvJ(c, 0, rhsPtr, nRhs, false);
00856 }
00857 }
00858 else
00859 {
00860 SUNDANCE_MSG2(dfVerb(), tab1 << "no derivatives to transform");
00861 }
00862 SUNDANCE_MSG2(dfVerb(), tab1 << "done transformations");
00863 }
00864
00865 jFlops = CellJacobianBatch::totalFlops() - jFlops;
00866 addFlops(flops + jFlops);
00867
00868 SUNDANCE_MSG2(dfVerb(),
00869 tab0 << "done QuadratureEvalMediator::fillFunctionCache()");
00870 }
00871
00872 void QuadratureEvalMediator::computePhysQuadPts() const
00873 {
00874 if (cacheIsValid())
00875 {
00876 Tabs tab0(0);
00877 SUNDANCE_MSG2(verb(), tab0 << "reusing phys quad points");
00878 }
00879 else
00880 {
00881 Tabs tab0(0);
00882 double jFlops = CellJacobianBatch::totalFlops();
00883 SUNDANCE_MSG2(verb(), tab0 << "computing phys quad points");
00884 physQuadPts_.resize(0);
00885 if (cellDim() != maxCellDim() && ElementIntegral::alwaysUseCofacets()
00886 && !isInternalBdry())
00887 {
00888 Tabs tab1;
00889 SUNDANCE_MSG2(verb(), tab1 << "using cofacets");
00890 SUNDANCE_MSG3(verb(), tab1 << "num cofacets = " << cofacetCellLID()->size());
00891 Array<Point> tmp;
00892 Array<int> cell(1);
00893 for (int c=0; c<cellLID()->size(); c++)
00894 {
00895 cell[0] = (*cofacetCellLID())[c];
00896 int facetIndex = facetIndices()[c];
00897 const Array<Point>& refFacetPts
00898 = (*(quadPtsReferredToMaxCell_.get(cellType())))[facetIndex];
00899 tmp.resize(refFacetPts.size());
00900 mesh().pushForward(maxCellDim(), cell, refFacetPts, tmp);
00901 for (int q=0; q<tmp.size(); q++)
00902 {
00903 physQuadPts_.append(tmp[q]);
00904 }
00905 }
00906 }
00907 else
00908 {
00909 const Array<Point>& refPts
00910 = *(quadPtsForReferenceCell_.get(cellType()));
00911 mesh().pushForward(cellDim(), *cellLID(),
00912 refPts, physQuadPts_);
00913 }
00914 addFlops(CellJacobianBatch::totalFlops() - jFlops);
00915 cacheIsValid() = true;
00916 SUNDANCE_OUT(this->verb() > 2,
00917 "phys quad: " << physQuadPts_);
00918 }
00919 }
00920
00921
00922 void QuadratureEvalMediator::print(std::ostream& os) const
00923 {
00924 if (cacheIsValid())
00925 {
00926 Tabs tab0;
00927 os << tab0 << "Physical quadrature points" << std::endl;
00928 int i=0;
00929 for (int c=0; c<cellLID()->size(); c++)
00930 {
00931 Tabs tab1;
00932 os << tab1 << "cell " << c << std::endl;
00933 for (int q=0; q<physQuadPts_.size()/cellLID()->size(); q++, i++)
00934 {
00935 Tabs tab2;
00936 os << tab2 << "q=" << q << " " << physQuadPts_[i] << std::endl;
00937 }
00938 }
00939 }
00940 }
00941
00942
00943 void QuadratureEvalMediator
00944 ::showResults(std::ostream& os,
00945 const RCP<SparsitySuperset>& sp,
00946 const Array<RCP<EvalVector> >& vecResults,
00947 const Array<double>& constantResults) const
00948 {
00949 Tabs tabs(0);
00950
00951
00952
00953
00954
00955 int maxlen = 25;
00956 for (int i=0; i<sp->numDerivs(); i++)
00957 {
00958 int s = sp->deriv(i).toString().length();
00959 if (s > maxlen) maxlen = s;
00960 }
00961
00962
00963 int vecIndex=0;
00964 int constIndex = 0;
00965 os << tabs << "Results Superset" << std::endl;
00966 for (int i=0; i<sp->numDerivs(); i++)
00967 {
00968 Tabs tab1;
00969 os << tab1 << i << " " ;
00970 os.width(maxlen);
00971 os.setf(std::ios_base::left, std::ios_base::adjustfield);
00972 os << sp->deriv(i).toString() ;
00973 Tabs tab2;
00974 switch(sp->state(i))
00975 {
00976 case ZeroDeriv:
00977 os << "Zero" << std::endl;
00978 break;
00979 case ConstantDeriv:
00980 os << "const val=" << constantResults[constIndex++] << std::endl;
00981 break;
00982 case VectorDeriv:
00983 if (vecResults[vecIndex].get()==0)
00984 {
00985 os << "{Null}";
00986 }
00987 else
00988 {
00989 const double* data = vecResults[vecIndex]->start();
00990 if (EvalVector::shadowOps())
00991 {
00992 TEUCHOS_TEST_FOR_EXCEPT(vecResults[vecIndex]->str().size()==0);
00993 os << vecResults[vecIndex]->str();
00994 }
00995 os << endl;
00996 int nQuad = numQuadPts(cellType());
00997 int k=0;
00998 TEUCHOS_TEST_FOR_EXCEPT(cellLID()->size() * nQuad
00999 != vecResults[vecIndex]->length());
01000 for (int c=0; c<cellLID()->size(); c++)
01001 {
01002 Tabs tab3;
01003 os << tab3 << "cell LID=" << (*cellLID())[c] << endl;
01004 for (int q=0; q<nQuad; q++, k++)
01005 {
01006 Tabs tab4;
01007 os << tab4 << "q=" << setw(5) << q
01008 << setw(10) << Utils::chop(data[k]) << endl;
01009 }
01010 }
01011 }
01012 vecIndex++;
01013 os << std::endl;
01014 break;
01015 }
01016 }
01017 }
