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 "SundanceCurveEvalMediator.hpp"
00032 #include "SundanceCoordExpr.hpp"
00033 #include "SundanceTempStack.hpp"
00034 #include "SundanceCellDiameterExpr.hpp"
00035 #include "SundanceCurveNormExpr.hpp"
00036 #include "SundanceCellVectorExpr.hpp"
00037 #include "SundanceSpatialDerivSpecifier.hpp"
00038 #include "SundanceDiscreteFunction.hpp"
00039 #include "SundanceDiscreteFuncElement.hpp"
00040 #include "SundanceCellJacobianBatch.hpp"
00041 #include "SundanceOut.hpp"
00042 #include "PlayaTabs.hpp"
00043 #include "PlayaExceptions.hpp"
00044 #include "SundanceCurveIntegralCalc.hpp"
00045
00046 #include "Teuchos_BLAS.hpp"
00047
00048
00049 using namespace Sundance;
00050 using namespace Teuchos;
00051 using namespace Playa;
00052
00053 TEUCHOS_TIMER(coordEvalTimer, "Quad mediator: coord eval")
00054
00055 CurveEvalMediator::CurveEvalMediator(
00056 const Mesh& mesh, const ParametrizedCurve& paramcurve ,
00057 int cellDim, const QuadratureFamily& quad) : StdFwkEvalMediator(mesh, cellDim) ,
00058 numEvaluationCases_(-1),
00059 quad_(quad),
00060 numQuadPtsForMaxCell_(0) ,
00061 paramcurve_(paramcurve) {
00062
00063
00064 maxCellType_ = mesh.cellType(mesh.spatialDim() );
00065 curveCellType_ = mesh.cellType(mesh.spatialDim()-1);
00066
00067
00068 Array<Point> quadPoints;
00069 Array<double> quadWeights;
00070 quad.getPoints(curveCellType_ , quadPoints , quadWeights );
00071 numQuadPtsForMaxCell_ = quadWeights.size();
00072 }
00073
00074 Time& CurveEvalMediator::coordEvaluationTimer()
00075 {
00076 return coordEvalTimer();
00077 }
00078
00079 void CurveEvalMediator::setCellType(const CellType& cellType,
00080 const CellType& maxCellType,
00081 bool isInternBdry)
00082 {
00083 StdFwkEvalMediator::setCellType(cellType, maxCellType, isInternBdry);
00084
00085 maxCellType_ = maxCellType;
00086
00087 Tabs tab;
00088
00089 SUNDANCE_MSG2(verb(), tab << "CurveEvalMediator::setCellType: cellType="
00090 << cellType << " cellDim=" << cellDim() << " maxCellType=" << maxCellType);
00091 SUNDANCE_MSG2(verb(), tab << "integration spec: ="
00092 << integrationCellSpec());
00093 SUNDANCE_MSG2(verb(), tab << "forbid cofacet integrations ="
00094 << forbidCofacetIntegrations());
00095 if (isInternalBdry())
00096 {
00097 SUNDANCE_MSG2(verb(), tab << "working on internal boundary");
00098 }
00099
00100
00101 if (cellType != maxCellType)
00102 {
00103
00104 }
00105 else
00106 {
00107 Tabs tab1;
00108 SUNDANCE_MSG2(verb(), tab1 << "no need for facet cases; work on original cell");
00109 numEvaluationCases_ = 1;
00110 }
00111
00112 SUNDANCE_MSG2(verb(), tab << "num eval cases =" << numEvaluationCases_);
00113
00114 }
00115
00116 int CurveEvalMediator::numQuadPts(const CellType& ct) const
00117 {
00118 if (ct != maxCellType_){
00119
00120 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "CurveEvalMediator::numQuadPts , cellType must be maxCellType_:" << ct);
00121 }
00122 return numQuadPtsForMaxCell_;
00123 }
00124
00125 void CurveEvalMediator::evalCellDiameterExpr(const CellDiameterExpr* expr,
00126 RCP<EvalVector>& vec) const
00127 {
00128
00129 Tabs tabs;
00130 SUNDANCE_MSG2(verb(),tabs
00131 << "CurveEvalMediator evaluating cell diameter expr "
00132 << expr->toString());
00133
00134 int nQuad = numQuadPts(cellType());
00135 int nCells = cellLID()->size();
00136
00137 SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad);
00138 Array<double> diameters;
00139 mesh().getCellDiameters(cellDim(), *cellLID(), diameters);
00140
00141 vec->resize(nQuad*nCells);
00142 double * const xx = vec->start();
00143 int k=0;
00144 for (int c=0; c<nCells; c++)
00145 {
00146 double h = diameters[c];
00147 for (int q=0; q<nQuad; q++, k++)
00148 {
00149 xx[k] = h;
00150 }
00151 }
00152 }
00153
00154 void CurveEvalMediator::evalCellVectorExpr(const CellVectorExpr* expr,
00155 RCP<EvalVector>& vec) const
00156 {
00157 Tabs tabs;
00158 SUNDANCE_MSG2(verb(),tabs
00159 << "CurveEvalMediator evaluating cell vector expr "
00160 << expr->toString());
00161
00162 int nQuad = numQuadPts(cellType());
00163 int nCells = cellLID()->size();
00164
00165 vec->resize(nQuad*nCells);
00166
00167 SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad);
00168 int dir = expr->componentIndex();
00169
00170 Array<Point> vectors;
00171 if (expr->isNormal())
00172 {
00173 mesh().outwardNormals(*cellLID(), vectors);
00174 }
00175 else
00176 {
00177 TEUCHOS_TEST_FOR_EXCEPTION(cellDim() != 1, std::runtime_error,
00178 "unable to compute tangent vectors for cell dim = " << cellDim());
00179 mesh().tangentsToEdges(*cellLID(), vectors);
00180 }
00181
00182 double * const xx = vec->start();
00183 int k=0;
00184 for (int c=0; c<nCells; c++)
00185 {
00186 double n = vectors[c][dir];
00187 for (int q=0; q<nQuad; q++, k++)
00188 {
00189 xx[k] = n;
00190 }
00191 }
00192 }
00193
00194 void CurveEvalMediator::evalCoordExpr(const CoordExpr* expr,
00195 RCP<EvalVector>& vec) const
00196 {
00197 Tabs tabs;
00198 SUNDANCE_MSG2(verb(),tabs
00199 << "CurveEvalMediator evaluating coord expr "
00200 << expr->toString());
00201
00202 TimeMonitor timer(coordEvalTimer());
00203
00204 int nQuad = numQuadPts(cellType());
00205 int nCells = cellLID()->size();
00206 int d = expr->dir();
00207
00208 SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad << ", nCells:" << nCells);
00209
00210
00211 vec->resize(nQuad*nCells);
00212 double * const xx = vec->start();
00213 int k=0;
00214 int maxDim = mesh().spatialDim() ;
00215 Array<int> cellLIDs_tmp(1);
00216 Array<Point> phyPoints;
00217 Array<Point> refPoints;
00218 Array<Point> refDevs;
00219 Array<Point> refNormal;
00220
00221
00222 for (int c=0; c<nCells; c++)
00223 {
00224 int maxCellLID = (*cellLID())[c];
00225 cellLIDs_tmp[0] = (*cellLID())[c];
00226
00227 if ( mesh().hasCurvePoints( maxCellLID , paramcurve_.myID() ))
00228 {
00229 mesh().getCurvePoints( maxCellLID , paramcurve_.myID() , refPoints , refDevs , refNormal );
00230 }
00231 else
00232 {
00233
00234 CurveIntegralCalc::getCurveQuadPoints( maxCellType_ , maxCellLID , mesh() , paramcurve_ , quad_ ,
00235 refPoints, refDevs , refNormal);
00236
00237 mesh().setCurvePoints( maxCellLID , paramcurve_.myID() , refPoints, refDevs , refNormal );
00238 }
00239
00240
00241 mesh().pushForward( maxDim , cellLIDs_tmp , refPoints , phyPoints );
00242
00243
00244 for (int q=0; q<nQuad; q++, k++)
00245 {
00246 xx[k] = phyPoints[q][d];
00247 }
00248 }
00249 }
00250
00251
00252 void CurveEvalMediator::evalCurveNormExpr(const CurveNormExpr* expr,
00253 RCP<EvalVector>& vec) const {
00254
00255 Tabs tabs;
00256 SUNDANCE_MSG2(verb(),tabs
00257 << "CurveEvalMediator evaluating curve normal expr "
00258 << expr->toString());
00259
00260 int nQuad = numQuadPts(cellType());
00261 int nCells = cellLID()->size();
00262 int d = expr->dir();
00263
00264 SUNDANCE_MSG3(verb(),tabs << "number of quad pts=" << nQuad);
00265
00266
00267 vec->resize(nCells*nQuad);
00268
00269 double * const xx = vec->start();
00270 int k=0;
00271 Array<int> cellLIDs_tmp(1);
00272 Array<Point> phyPoints;
00273 Array<Point> refPoints;
00274 Array<Point> refDevs;
00275 Array<Point> refNormal;
00276
00277
00278 for (int c=0; c<nCells; c++)
00279 {
00280 int maxCellLID = (*cellLID())[c];
00281 cellLIDs_tmp[0] = (*cellLID())[c];
00282
00283 if ( mesh().hasCurvePoints( maxCellLID , paramcurve_.myID() ))
00284 {
00285 mesh().getCurvePoints( maxCellLID , paramcurve_.myID() , refPoints , refDevs , refNormal );
00286 }
00287 else
00288 {
00289
00290 CurveIntegralCalc::getCurveQuadPoints( maxCellType_ , maxCellLID , mesh() , paramcurve_ , quad_ ,
00291 refPoints, refDevs , refNormal);
00292
00293 mesh().setCurvePoints( maxCellLID , paramcurve_.myID() , refPoints, refDevs , refNormal );
00294 }
00295
00296
00297 for (int q=0; q<nQuad; q++, k++)
00298 {
00299 xx[k] = refNormal[q][d];
00300 }
00301 }
00302 }
00303
00304
00305 void CurveEvalMediator
00306 ::evalDiscreteFuncElement(const DiscreteFuncElement* expr,
00307 const Array<MultiIndex>& multiIndices,
00308 Array<RCP<EvalVector> >& vec) const
00309 {
00310
00311 int verbo = dfVerb();
00312 Tabs tab1;
00313 SUNDANCE_MSG2(verbo , tab1
00314 << "CurveEvalMediator evaluating Discrete Function expr " << expr->toString());
00315
00316 const DiscreteFunctionData* f = DiscreteFunctionData::getData(expr);
00317
00318 TEUCHOS_TEST_FOR_EXCEPTION(f==0, std::logic_error,
00319 "QuadratureEvalMediator::evalDiscreteFuncElement() called "
00320 "with expr that is not a discrete function");
00321
00322 SUNDANCE_MSG2(verbo , tab1 << "After casting DiscreteFunctionData" << expr->toString());
00323
00324 RCP<Array<Array<double> > > localValues;
00325 Array<int> cellLIDs_tmp(1);
00326 Array<Point> phyPoints;
00327 Array<Point> refPoints;
00328 Array<Point> refDevs;
00329 Array<Point> refNormal;
00330 int nCells = cellLID()->size();
00331
00332 RCP<const MapStructure> mapStruct;
00333 int myIndex = expr->myIndex();
00334 int nQuad = numQuadPtsForMaxCell_;
00335 Array<int> k(multiIndices.size(),0);
00336
00337 Teuchos::BLAS<int,double> blas;
00338
00339 SUNDANCE_MSG2(verbo , tab1 << "After declaring BLAS: " << expr->toString());
00340
00341
00342 for (int i=0; i<multiIndices.size(); i++)
00343 {
00344 vec[i]->resize(nCells*nQuad);
00345 }
00346
00347
00348 for (int c=0; c<nCells; c++)
00349 {
00350 int maxCellLID = (*cellLID())[c];
00351 localValues = rcp(new Array<Array<double> >());
00352 SUNDANCE_MSG2(verbo , tab1 << "Cell:" << c << " of " << nCells << " , maxCellLID:" << maxCellLID );
00353
00354 cellLIDs_tmp[0] = (*cellLID())[c];
00355 SUNDANCE_MSG2(verbo , tab1 << " Before calling f->getLocalValues:" << cellLIDs_tmp.size() <<
00356 " f==0 : " << (f==0) << " tmp:" << f->mesh().spatialDim());
00357
00358 mapStruct = f->getLocalValues(maxCellDim(), cellLIDs_tmp , *localValues);
00359 SUNDANCE_MSG2(verbo , tab1 << " After getting mapStruct:" << maxCellLID );
00360 SUNDANCE_MSG2(verbo , tab1 << " mapStruct->numBasisChunks():" << mapStruct->numBasisChunks() );
00361 int chunk = mapStruct->chunkForFuncID(myIndex);
00362 int funcIndex = mapStruct->indexForFuncID(myIndex);
00363 int nFuncs = mapStruct->numFuncs(chunk);
00364 SUNDANCE_MSG2(verbo , tab1 << " chunk:" << chunk );
00365 SUNDANCE_MSG2(verbo , tab1 << " funcIndex:" << funcIndex );
00366 SUNDANCE_MSG2(verbo , tab1 << " nFuncs:" << nFuncs );
00367
00368
00369 BasisFamily basis = rcp_dynamic_cast<BasisFamilyBase>(mapStruct->basis(chunk));
00370 int nNodesTotal = basis.nReferenceDOFsWithFacets(maxCellType(), maxCellType());
00371
00372
00373 if ( mesh().hasCurvePoints( maxCellLID , paramcurve_.myID() ))
00374 {
00375 mesh().getCurvePoints( maxCellLID , paramcurve_.myID() , refPoints , refDevs , refNormal );
00376 }
00377 else
00378 {
00379
00380 CurveIntegralCalc::getCurveQuadPoints( maxCellType_ , maxCellLID , mesh() , paramcurve_ , quad_ ,
00381 refPoints, refDevs , refNormal);
00382
00383 mesh().setCurvePoints( maxCellLID , paramcurve_.myID() , refPoints, refDevs , refNormal );
00384 }
00385
00386
00387 SUNDANCE_MSG2(verbo , tab1 << " multiIndices.size()" << multiIndices.size() );
00388 for (int i=0; i<multiIndices.size(); i++)
00389 {
00390
00391 int nDerivResults = 1;
00392 if ( multiIndices[i].order() == 1 ) nDerivResults = maxCellDim();
00393
00394 int pDir = 0;
00395 int derivNum = 1;
00396 MultiIndex mi;
00397 SUNDANCE_MSG2(verbo , tab1 << " before asking anything i = " << i);
00398 SUNDANCE_MSG2(verbo , tab1 << " multiindex order : " << multiIndices[i].order());
00399 SUNDANCE_MSG2(verbo , tab1 << " multiindex : " << multiIndices[i] );
00400 if (multiIndices[i].order() > 0){
00401 pDir = multiIndices[i].firstOrderDirection();
00402 mi[pDir] = 1;
00403 derivNum = mesh().spatialDim();
00404 }
00405
00406 Array<Array<double> > result(nQuad*derivNum);
00407 Array<Array<Array<double> > > tmp;
00408
00409 int offs = nNodesTotal;
00410
00411
00412 for (int deriv = 0 ; deriv < derivNum ; deriv++)
00413 {
00414
00415 if (multiIndices[i].order() > 0){
00416
00417 MultiIndex mi_tmp;
00418 mi_tmp[deriv] = 1;
00419 SpatialDerivSpecifier deriv(mi_tmp);
00420 SUNDANCE_MSG2(verbo , tab1 << "computing derivatives : " << deriv << " on reference cell ");
00421 basis.refEval( maxCellType_ , refPoints , deriv, tmp , verbo );
00422 } else {
00423 SpatialDerivSpecifier deriv(mi);
00424
00425 SUNDANCE_MSG2(verbo , tab1 << "computing values reference cell ");
00426 basis.refEval( maxCellType_ , refPoints , deriv, tmp , verbo );
00427 }
00428
00429 SUNDANCE_MSG2(verbo , tab1 << "resize result vector , offs:" << offs);
00430 for (int q=0; q<nQuad; q++){
00431 result[nQuad*deriv + q].resize(offs);
00432 }
00433
00434
00435 SUNDANCE_MSG2(verbo , tab1 << "copy results ");
00436 int offs1 = 0;
00437 for (int q=0; q<nQuad; q++)
00438 {
00439 offs1 = 0;
00440 for (int d=0; d<basis.dim(); d++)
00441 {
00442 int nNodes = tmp[d][q].size();
00443 for (int n=0; n<nNodes; n++ , offs1++ )
00444 {
00445 result[nQuad*deriv + q][offs1] = tmp[d][q][n];
00446 }
00447 }
00448 }
00449 }
00450
00451
00452 SUNDANCE_MSG2(verbo , tab1 << "summing up values , funcIndex:" << funcIndex << " offs:" << offs);
00453 for (int deriv = 0 ; deriv < derivNum ; deriv++)
00454 {
00455 for (int q=0; q<nQuad; q++)
00456 {
00457 double sum = 0.0;
00458
00459 for (int n = 0 ; n < offs ; n++){
00460 sum = sum + result[nQuad*deriv + q][n] * (*localValues)[chunk][funcIndex*offs + n];
00461 }
00462
00463 result[nQuad*deriv + q][0] = sum;
00464 }
00465 }
00466
00467
00468 const CellJacobianBatch& J = JTrans();
00469 if (mi.order()==1)
00470 {
00471 Tabs tab1;
00472 Tabs tab2;
00473 SUNDANCE_MSG2(verbo, tab2 << "Jacobian batch nCells=" << J.numCells());
00474 SUNDANCE_MSG2(verbo, tab2 << "Jacobian batch cell dim=" << J.cellDim());
00475 SUNDANCE_MSG2(verbo, tab2 << "Jacobian batch spatial dim=" << J.spatialDim());
00476
00477 Array<double> invJ;
00478 J.getInvJ(c, invJ);
00479 for (int q=0; q<nQuad; q++)
00480 {
00481 double sum = 0.0;
00482 for (int deriv = 0 ; deriv < derivNum ; deriv++)
00483 {
00484
00485 sum = sum + result[nQuad*deriv + q][0] * invJ[derivNum*pDir + deriv];
00486 }
00487
00488 result[q][0] = sum;
00489 }
00490 }
00491
00492
00493
00494 double* vecPtr = vec[i]->start();
00495 for (int q=0; q<nQuad; q++, k[i]++)
00496 {
00497 vecPtr[k[i]] = result[q][0];
00498 }
00499 SUNDANCE_MSG2(verbo , tab1 << " END copy results back ");
00500 }
00501 SUNDANCE_MSG2(verbo , tab1 << " END loop over multiindex ");
00502 }
00503 SUNDANCE_MSG2(verbo , tab1 << " END loop over cells ");
00504 }
00505
00506
00507 void CurveEvalMediator::print(std::ostream& os) const
00508 {
00509
00510 }
00511
00512