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 "SundanceAssembler.hpp"
00032 #include "SundanceDOFMapBuilder.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "PlayaTabs.hpp"
00035 #include "SundanceCellFilter.hpp"
00036 #include "SundanceCellSet.hpp"
00037 #include "SundanceTrivialGrouper.hpp"
00038 #include "SundanceDOFMapBase.hpp"
00039 #include "SundanceEquationSet.hpp"
00040 #include "SundanceDiscreteSpace.hpp"
00041 #include "SundanceDiscreteFunction.hpp"
00042 #include "SundanceIntegralGroup.hpp"
00043 #include "SundanceGrouperBase.hpp"
00044 #include "SundanceEvalManager.hpp"
00045 #include "SundanceStdFwkEvalMediator.hpp"
00046 #include "SundanceEvaluatableExpr.hpp"
00047 #include "PlayaLoadableVector.hpp"
00048 #include "PlayaLoadableMatrix.hpp"
00049 #include "SundanceQuadratureEvalMediator.hpp"
00050 #include "SundanceCurveEvalMediator.hpp"
00051 #include "SundanceEvaluator.hpp"
00052 #include "Teuchos_Time.hpp"
00053 #include "Teuchos_TimeMonitor.hpp"
00054 #include "Epetra_HashTable.h"
00055 #include "SundanceIntHashSet.hpp"
00056 #include "PlayaDefaultBlockVectorSpaceDecl.hpp"
00057 #include "PlayaBlockOperatorBaseDecl.hpp"
00058 #include "PlayaSimpleBlockOpDecl.hpp"
00059 #include "SundanceAssemblyKernelBase.hpp"
00060 #include "SundanceVectorAssemblyKernel.hpp"
00061 #include "SundanceMatrixVectorAssemblyKernel.hpp"
00062 #include "SundanceFunctionalAssemblyKernel.hpp"
00063 #include "SundanceFunctionalGradientAssemblyKernel.hpp"
00064 #include "SundanceAssemblyTransformationBuilder.hpp"
00065 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00066 #include "PlayaLinearOperatorImpl.hpp"
00067 #include "PlayaSimpleBlockOpImpl.hpp"
00068 #include "PlayaDefaultBlockVectorImpl.hpp"
00069 #endif
00070
00071
00072
00073 using namespace Teuchos;
00074 using namespace Playa;
00075 using std::max;
00076 using std::min;
00077
00078
00079 static Time& assemblerCtorTimer()
00080 {
00081 static RCP<Time> rtn
00082 = TimeMonitor::getNewTimer("assembler ctor");
00083 return *rtn;
00084 }
00085
00086
00087
00088
00089 static Time& graphBuildTimer()
00090 {
00091 static RCP<Time> rtn
00092 = TimeMonitor::getNewTimer("matrix graph determination");
00093 return *rtn;
00094 }
00095
00096 static Time& colSearchTimer()
00097 {
00098 static RCP<Time> rtn
00099 = TimeMonitor::getNewTimer("graph column processing");
00100 return *rtn;
00101 }
00102
00103
00104
00105 static Time& matAllocTimer()
00106 {
00107 static RCP<Time> rtn
00108 = TimeMonitor::getNewTimer("matrix allocation");
00109 return *rtn;
00110 }
00111
00112 static Time& matFinalizeTimer()
00113 {
00114 static RCP<Time> rtn
00115 = TimeMonitor::getNewTimer("matrix graph packing");
00116 return *rtn;
00117 }
00118
00119 static Time& graphFlatteningTimer()
00120 {
00121 static RCP<Time> rtn
00122 = TimeMonitor::getNewTimer("tmp graph flattening");
00123 return *rtn;
00124 }
00125
00126
00127
00128 Assembler
00129 ::Assembler(const Mesh& mesh,
00130 const RCP<EquationSet>& eqn,
00131 const Array<VectorType<double> >& rowVectorType,
00132 const Array<VectorType<double> >& colVectorType,
00133 bool partitionBCs)
00134 : partitionBCs_(partitionBCs),
00135 matNeedsConfiguration_(true),
00136 matNeedsFinalization_(true),
00137 numConfiguredColumns_(0),
00138 mesh_(mesh),
00139 eqn_(eqn),
00140 rowMap_(),
00141 colMap_(),
00142 externalRowSpace_(eqn->numVarBlocks()),
00143 externalColSpace_(eqn->numUnkBlocks()),
00144 privateRowSpace_(eqn->numVarBlocks()),
00145 privateColSpace_(eqn->numUnkBlocks()),
00146 bcRows_(eqn->numVarBlocks()),
00147 bcCols_(eqn->numUnkBlocks()),
00148 rqc_(),
00149 contexts_(),
00150 isBCRqc_(),
00151 isInternalBdry_(),
00152 groups_(),
00153 mediators_(),
00154 evalExprs_(),
00155 evalMgr_(rcp(new EvalManager())),
00156 isBCRow_(eqn->numVarBlocks()),
00157 isBCCol_(eqn->numUnkBlocks()),
00158 remoteBCCols_(eqn->numUnkBlocks()),
00159 lowestRow_(eqn->numVarBlocks()),
00160 lowestCol_(eqn->numUnkBlocks()),
00161 rowVecType_(rowVectorType),
00162 colVecType_(colVectorType),
00163 testIDToBlockMap_(),
00164 unkIDToBlockMap_()
00165 {
00166 TimeMonitor timer(assemblerCtorTimer());
00167 init(mesh, eqn);
00168 }
00169
00170 Assembler
00171 ::Assembler(const Mesh& mesh,
00172 const RCP<EquationSet>& eqn)
00173 : partitionBCs_(false),
00174 matNeedsConfiguration_(true),
00175 matNeedsFinalization_(true),
00176 numConfiguredColumns_(0),
00177 mesh_(mesh),
00178 eqn_(eqn),
00179 rowMap_(),
00180 colMap_(),
00181 externalRowSpace_(eqn->numVarBlocks()),
00182 externalColSpace_(eqn->numUnkBlocks()),
00183 privateRowSpace_(eqn->numVarBlocks()),
00184 privateColSpace_(eqn->numUnkBlocks()),
00185 bcRows_(eqn->numVarBlocks()),
00186 bcCols_(eqn->numUnkBlocks()),
00187 rqc_(),
00188 contexts_(),
00189 isBCRqc_(),
00190 isInternalBdry_(),
00191 groups_(),
00192 mediators_(),
00193 evalExprs_(),
00194 evalMgr_(rcp(new EvalManager())),
00195 isBCRow_(eqn->numVarBlocks()),
00196 isBCCol_(eqn->numUnkBlocks()),
00197 remoteBCCols_(eqn->numUnkBlocks()),
00198 lowestRow_(eqn->numVarBlocks()),
00199 lowestCol_(eqn->numUnkBlocks()),
00200 rowVecType_(),
00201 colVecType_(),
00202 testIDToBlockMap_(),
00203 unkIDToBlockMap_(),
00204 fixedParamIDToVectorNumber_()
00205 {
00206 TimeMonitor timer(assemblerCtorTimer());
00207 init(mesh, eqn);
00208 }
00209
00210 void Assembler::init(const Mesh& mesh, const RCP<EquationSet>& eqn)
00211 {
00212 Tabs tab0(0);
00213
00214
00215 int verb = 0;
00216 if (eqn->hasActiveWatchFlag())
00217 verb = eqn->maxWatchFlagSetting("assembler setup");
00218
00219 SUNDANCE_BANNER1(verb, tab0, "Assembler setup");
00220
00221
00222
00223 RCP<GrouperBase> grouper
00224 = rcp(new TrivialGrouper());
00225
00226
00227
00228
00229 const Set<ComputationType>& compTypes = eqn->computationTypes();
00230
00231
00232
00233 DOFMapBuilder mapBuilder(eqn->maxWatchFlagSetting("dof map setup"));
00234
00235 if (compTypes.contains(VectorOnly)
00236 || compTypes.contains(Sensitivities)
00237 || compTypes.contains(FunctionalAndGradient))
00238 {
00239 Tabs tab1;
00240 SUNDANCE_MSG2(verb, tab1 << "building row spaces");
00241 mapBuilder = DOFMapBuilder(mesh, eqn->fsr(), partitionBCs_,
00242 eqn->maxWatchFlagSetting("dof map setup"));
00243
00244 rowMap_ = mapBuilder.rowMap();
00245 isBCRow_ = mapBuilder.isBCRow();
00246 isBCCol_ = mapBuilder.isBCCol();
00247 lowestRow_.resize(eqn_->numVarBlocks());
00248
00249 for (int b=0; b<eqn_->numVarBlocks(); b++)
00250 {
00251 Tabs tab2;
00252 lowestRow_[b] = rowMap_[b]->lowestLocalDOF();
00253 SUNDANCE_MSG2(verb, tab2 << "block " << b << ": lowest row="
00254 << lowestRow_[b]);
00255 externalRowSpace_[b] = rcp(
00256 new DiscreteSpace(mesh, testBasisArray(mapBuilder.fsr())[b],
00257 rowMap_[b], rowVecType_[b]));
00258 privateRowSpace_[b] = externalRowSpace_[b];
00259 SUNDANCE_MSG2(verb, tab2 << "block " << b << ": done forming row space");
00260 }
00261 }
00262
00263 if (!eqn->isFunctionalCalculator())
00264 {
00265 Tabs tab1;
00266
00267 SUNDANCE_MSG2(verb, tab1 << "building column spaces");
00268 colMap_ = mapBuilder.colMap();
00269
00270 for (int b=0; b<eqn_->numUnkBlocks(); b++)
00271 {
00272 Tabs tab2;
00273 externalColSpace_[b]
00274 = rcp(new DiscreteSpace(mesh, unkBasisArray(mapBuilder.fsr())[b],
00275 colMap_[b], colVecType_[b]));
00276 privateColSpace_[b] = externalColSpace_[b];
00277 SUNDANCE_MSG2(verb, tab2 << "block " << b << ": done forming col space");
00278 }
00279
00280
00281
00282 groups_.put(MatrixAndVector, Array<Array<RCP<IntegralGroup> > >());
00283 rqcRequiresMaximalCofacets_.put(MatrixAndVector,
00284 Array<IntegrationCellSpecifier>());
00285 skipRqc_.put(MatrixAndVector, Array<int>());
00286 contexts_.put(MatrixAndVector, Array<EvalContext>());
00287 evalExprs_.put(MatrixAndVector, Array<const EvaluatableExpr*>());
00288
00289
00290 groups_.put(VectorOnly, Array<Array<RCP<IntegralGroup> > >());
00291 rqcRequiresMaximalCofacets_.put(VectorOnly,
00292 Array<IntegrationCellSpecifier>());
00293 skipRqc_.put(VectorOnly, Array<int>());
00294 contexts_.put(VectorOnly, Array<EvalContext>());
00295 evalExprs_.put(VectorOnly, Array<const EvaluatableExpr*>());
00296
00297 if (eqn->isSensitivityCalculator())
00298 {
00299 fixedParamIDToVectorNumber_
00300 = eqn->fsr()->fixedParamIDToReducedFixedParamIDMap();
00301
00302
00303 groups_.put(Sensitivities, Array<Array<RCP<IntegralGroup> > >());
00304 rqcRequiresMaximalCofacets_.put(Sensitivities,
00305 Array<IntegrationCellSpecifier>());
00306 skipRqc_.put(Sensitivities, Array<int>());
00307 contexts_.put(Sensitivities, Array<EvalContext>());
00308 evalExprs_.put(Sensitivities, Array<const EvaluatableExpr*>());
00309
00310 }
00311 }
00312 else
00313 {
00314
00315 groups_.put(FunctionalAndGradient, Array<Array<RCP<IntegralGroup> > >());
00316 rqcRequiresMaximalCofacets_.put(FunctionalAndGradient,
00317 Array<IntegrationCellSpecifier>());
00318 skipRqc_.put(FunctionalAndGradient, Array<int>());
00319 contexts_.put(FunctionalAndGradient, Array<EvalContext>());
00320 evalExprs_.put(FunctionalAndGradient, Array<const EvaluatableExpr*>());
00321
00322 groups_.put(FunctionalOnly, Array<Array<RCP<IntegralGroup> > >());
00323 rqcRequiresMaximalCofacets_.put(FunctionalOnly,
00324 Array<IntegrationCellSpecifier>());
00325 skipRqc_.put(FunctionalOnly, Array<int>());
00326 contexts_.put(FunctionalOnly, Array<EvalContext>());
00327 evalExprs_.put(FunctionalOnly, Array<const EvaluatableExpr*>());
00328 }
00329
00330
00331
00332
00333
00334
00335 SUNDANCE_MSG1(verb, tab0 << std::endl
00336 << tab0 << "=== setting up non-BC region-quadrature combinations");
00337
00338 for (int r=0; r<eqn->regionQuadCombos().size(); r++)
00339 {
00340 Tabs tab1;
00341 Tabs tab12;
00342 const RegionQuadCombo& rqc = eqn->regionQuadCombos()[r];
00343
00344
00345 bool watchMe = rqc.watch().isActive();
00346
00347 int rqcVerb = verb;
00348 int integralCtorVerb = 0;
00349 int integrationVerb = 0;
00350 int integralTransformVerb = 0;
00351 if (watchMe)
00352 {
00353 rqcVerb = max(4,rqcVerb);
00354 integralCtorVerb = rqc.watch().param("integration setup");
00355 integrationVerb = rqc.watch().param("integration");
00356 integralTransformVerb = rqc.watch().param("integral transformation");
00357 }
00358
00359
00360
00361 rqc_.append(rqc);
00362 isBCRqc_.append(false);
00363
00364
00365 const Expr& expr = eqn->expr(rqc);
00366
00367 SUNDANCE_MSG2(rqcVerb, tab1 << std::endl << tab1 << "------------------------------------------------");
00368 SUNDANCE_MSG2(rqcVerb, tab1 << "initializing assembly for"
00369 << std::endl << tab12 << "r=" << r << " rqc="
00370 << rqc << std::endl << tab12 << std::endl << tab12 << "------------------------------"
00371 << std::endl << tab12 << "expr = " << expr
00372 << std::endl << tab12 << "------------------------------"
00373 );
00374
00375
00376
00377 int cellDim = CellFilter(rqc.domain()).dimension(mesh);
00378 CellType cellType = mesh.cellType(cellDim);
00379 CellType maxCellType = mesh.cellType(mesh.spatialDim());
00380 QuadratureFamily quad(rqc.quad());
00381
00382
00383 bool isInternalBdry = detectInternalBdry(cellDim, rqc.domain());
00384 isInternalBdry_.append(isInternalBdry);
00385
00386 SUNDANCE_MSG2(rqcVerb, tab12 << "isInternalBdry=" << isInternalBdry);
00387
00388
00389 bool rqcUsed = false;
00390
00391 for (Set<ComputationType>::const_iterator
00392 i=eqn->computationTypes().begin();
00393 i!=eqn->computationTypes().end();
00394 i++)
00395 {
00396 Tabs tab2;
00397 const ComputationType& compType = *i;
00398 SUNDANCE_MSG2(rqcVerb, tab12 << std::endl << tab12
00399 << "** computation type " << compType);
00400
00401
00402
00403
00404 if (eqn->skipRqc(compType, rqc))
00405 {
00406 SUNDANCE_MSG2(rqcVerb, tab2 << "RQC not needed for computation type "
00407 << compType);
00408 skipRqc_[compType].append(true);
00409 EvalContext dummy;
00410 const EvaluatableExpr* dummyEx = 0;
00411 Array<RCP<IntegralGroup> > dummyGroups;
00412 IntegrationCellSpecifier dummyCellSpec ;
00413 contexts_[compType].append(dummy);
00414 evalExprs_[compType].append(dummyEx);
00415 groups_[compType].append(dummyGroups);
00416 rqcRequiresMaximalCofacets_[compType].append(dummyCellSpec);
00417 }
00418 else
00419 {
00420
00421
00422 rqcUsed = true;
00423 skipRqc_[compType].append(false);
00424
00425
00426
00427 EvalContext context = eqn->rqcToContext(compType, rqc);
00428
00429 SUNDANCE_MSG2(rqcVerb, tab2 << "context " << context.brief());
00430
00431
00432 contexts_[compType].append(context);
00433
00434
00435 const EvaluatableExpr* ee = EvaluatableExpr::getEvalExpr(expr);
00436 evalExprs_[compType].append(ee);
00437
00438
00439
00440
00441 const RCP<SparsitySuperset>& sparsity
00442 = ee->sparsitySuperset(context);
00443 SUNDANCE_MSG3(rqcVerb, tab2 << "sparsity pattern " << *sparsity);
00444
00445
00446
00447 Array<RCP<IntegralGroup> > groups;
00448 grouper->setVerb(integralCtorVerb, integrationVerb, integralTransformVerb);
00449 grouper->findGroups(*eqn, maxCellType, mesh.spatialDim(),
00450 cellType, cellDim, quad, sparsity, isInternalBdry, groups , rqc.paramCurve() , mesh_ );
00451 grouper->setVerb(0,0,0);
00452 groups_[compType].append(groups);
00453
00454
00455
00456 IntegrationCellSpecifier cellSpec
00457 = whetherToUseCofacets(groups, ee,
00458 cellDim==mesh_.spatialDim(), rqcVerb);
00459 SUNDANCE_MSG2(rqcVerb, tab2 << "integration: " << cellSpec);
00460 rqcRequiresMaximalCofacets_[compType].append(cellSpec);
00461 }
00462 SUNDANCE_MSG2(rqcVerb, tab12
00463 << "done with computation type " << compType);
00464 }
00465
00466
00467
00468
00469
00470
00471 if (rqcUsed)
00472 {
00473 SUNDANCE_MSG2(rqcVerb, tab12 << "creating evaluation mediator for rqc="
00474 << rqc);
00475 SUNDANCE_MSG2(rqcVerb, tab12 << "expr = " << expr);
00476
00477
00478
00479
00480
00481
00482 if ( rqc.paramCurve().isCurveIntegral() && rqc.paramCurve().isCurveValid() )
00483 {
00484 mediators_.append(rcp(new CurveEvalMediator(mesh, rqc.paramCurve() , cellDim, quad)));
00485 }
00486 else
00487 {
00488 mediators_.append(rcp(new QuadratureEvalMediator(mesh, cellDim, quad)));
00489 }
00490 }
00491 else
00492 {
00493 SUNDANCE_MSG2(rqcVerb, tab1 << "creating empty eval mediator for unused rqc");
00494 mediators_.append(RCP<StdFwkEvalMediator>());
00495 }
00496 SUNDANCE_MSG2(rqcVerb, tab1
00497 << "done with RQC");
00498 }
00499
00500
00501
00502
00503 SUNDANCE_MSG1(verb, tab0 << std::endl
00504 << tab0 << "=== setting up BC region-quadrature combinations");
00505
00506 for (int r=0; r<eqn->bcRegionQuadCombos().size(); r++)
00507 {
00508 Tabs tab1;
00509 const RegionQuadCombo& rqc = eqn->bcRegionQuadCombos()[r];
00510
00511
00512 bool watchMe = rqc.watch().isActive();
00513 int rqcVerb = verb;
00514 int integralCtorVerb = 0;
00515 int integrationVerb = 0;
00516 int integralTransformVerb = 0;
00517 if (watchMe)
00518 {
00519 rqcVerb = max(4,rqcVerb);
00520 integralCtorVerb = rqc.watch().param("integration setup");
00521 integrationVerb = rqc.watch().param("integration");
00522 integralTransformVerb = rqc.watch().param("integral transformation");
00523 }
00524
00525
00526 rqc_.append(rqc);
00527 isBCRqc_.append(true);
00528
00529
00530
00531 const Expr& expr = eqn->bcExpr(rqc);
00532
00533 SUNDANCE_MSG2(rqcVerb, tab1 << std::endl << tab1
00534 << "------------------------------------------------");
00535 SUNDANCE_MSG1(rqcVerb, tab1 << "initializing assembly for BC r=" << r
00536 << " rqc="
00537 << rqc << std::endl << tab1
00538 << "expr = " << expr);
00539
00540
00541 int cellDim = CellFilter(rqc.domain()).dimension(mesh);
00542 CellType cellType = mesh.cellType(cellDim);
00543 CellType maxCellType = mesh.cellType(mesh.spatialDim());
00544 QuadratureFamily quad(rqc.quad());
00545
00546
00547 bool isInternalBdry = detectInternalBdry(cellDim, rqc.domain());
00548 isInternalBdry_.append(isInternalBdry);
00549
00550
00551 bool rqcUsed = false;
00552
00553 for (Set<ComputationType>::const_iterator
00554 i=eqn->computationTypes().begin();
00555 i!=eqn->computationTypes().end();
00556 i++)
00557 {
00558 Tabs tab2;
00559 const ComputationType& compType = *i;
00560 SUNDANCE_MSG2(rqcVerb, tab1 << std::endl << tab1
00561 << "** computation type " << compType);
00562
00563
00564
00565
00566
00567 if (eqn->skipBCRqc(compType, rqc))
00568 {
00569 SUNDANCE_MSG2(rqcVerb,
00570 tab2 << "this rqc not needed for computation type " << compType);
00571 skipRqc_[compType].append(true);
00572 EvalContext dummy;
00573 const EvaluatableExpr* dummyEx = 0;
00574 Array<RCP<IntegralGroup> > dummyGroups;
00575 IntegrationCellSpecifier dummyCellSpec ;
00576 contexts_[compType].append(dummy);
00577 evalExprs_[compType].append(dummyEx);
00578 groups_[compType].append(dummyGroups);
00579 rqcRequiresMaximalCofacets_[compType].append(dummyCellSpec);
00580 }
00581 else
00582 {
00583
00584
00585 rqcUsed = true;
00586 skipRqc_[compType].append(false);
00587
00588
00589
00590 EvalContext context = eqn->bcRqcToContext(compType, rqc);
00591 SUNDANCE_MSG2(rqcVerb, tab2 << "context " << context);
00592
00593
00594 contexts_[compType].append(context);
00595 const EvaluatableExpr* ee = EvaluatableExpr::getEvalExpr(expr);
00596 evalExprs_[compType].append(ee);
00597 const RCP<SparsitySuperset>& sparsity
00598 = ee->sparsitySuperset(context);
00599 SUNDANCE_MSG3(rqcVerb, tab2 << "sparsity pattern " << *sparsity);
00600
00601 Array<RCP<IntegralGroup> > groups;
00602 grouper->setVerb(integralCtorVerb, integrationVerb, integralTransformVerb);
00603 grouper->findGroups(*eqn, maxCellType, mesh.spatialDim(),
00604 cellType, cellDim, quad, sparsity, isInternalBdry, groups , rqc.paramCurve() , mesh_ );
00605 grouper->setVerb(0,0,0);
00606 groups_[compType].append(groups);
00607 IntegrationCellSpecifier cellSpec
00608 = whetherToUseCofacets(groups, ee,
00609 cellDim==mesh_.spatialDim(), rqcVerb);
00610 SUNDANCE_MSG2(rqcVerb, tab2 << "integration: " << cellSpec);
00611 rqcRequiresMaximalCofacets_[compType].append(cellSpec);
00612 SUNDANCE_MSG2(rqcVerb, tab1
00613 << "done with computation type " << compType);
00614 }
00615
00616
00617
00618
00619 if (rqcUsed)
00620 {
00621 SUNDANCE_MSG2(rqcVerb, tab1 << "creating evaluation mediator for BC rqc="
00622 << rqc << std::endl << tab1 << "expr = " << expr);
00623
00624 if ( rqc.paramCurve().isCurveIntegral() && rqc.paramCurve().isCurveValid() )
00625 {
00626 mediators_.append(rcp(new CurveEvalMediator(mesh, rqc.paramCurve(), cellDim, quad)));
00627 }
00628 else
00629 {
00630 mediators_.append(rcp(new QuadratureEvalMediator(mesh, cellDim, quad)));
00631 }
00632 }
00633 else
00634 {
00635 SUNDANCE_MSG2(rqcVerb, tab1 << "creating empty eval mediator for unused BC rqc");
00636 mediators_.append(RCP<StdFwkEvalMediator>());
00637 }
00638 }
00639 SUNDANCE_MSG2(rqcVerb, tab1
00640 << "done with BC RQC");
00641 }
00642 }
00643
00644 bool Assembler::detectInternalBdry(int cellDim,
00645 const CellFilter& filter) const
00646 {
00647 int d = mesh_.spatialDim();
00648 if (cellDim == d-1)
00649 {
00650 CellSet cells = filter.getCells(mesh_);
00651 for (CellIterator c=cells.begin(); c!=cells.end(); c++)
00652 {
00653 if (mesh_.numMaxCofacets(cellDim, *c) > 1) return true;
00654 }
00655 }
00656 return false;
00657 }
00658
00659 IntegrationCellSpecifier Assembler::whetherToUseCofacets(
00660 const Array<RCP<IntegralGroup> >& groups,
00661 const EvaluatableExpr* ee,
00662 bool isMaximalCell,
00663 int verb) const
00664 {
00665 Tabs tab;
00666 SUNDANCE_MSG2(verb,
00667 tab << "deciding whether to use cofacet cells for some integrations");
00668
00669 if (isMaximalCell)
00670 {
00671 Tabs tab1;
00672 SUNDANCE_MSG2(verb,
00673 tab1 << "cofacets not needed because cells are maximal");
00674 return NoTermsNeedCofacets;
00675 }
00676
00677 IntegrationCellSpecifier cellSpec = SomeTermsNeedCofacets;
00678
00679 bool allTermsNeedCofacets = true;
00680 bool noTermsNeedCofacets = true;
00681 for (int g=0; g<groups.size(); g++)
00682 {
00683 Tabs tab1;
00684 switch(groups[g]->usesMaximalCofacets())
00685 {
00686 case NoTermsNeedCofacets:
00687 allTermsNeedCofacets = false;
00688 SUNDANCE_MSG2(verb,
00689 tab1 << "integral group " << g << " does not need cofacets");
00690 break;
00691 case AllTermsNeedCofacets:
00692 case SomeTermsNeedCofacets:
00693 noTermsNeedCofacets = false;
00694 SUNDANCE_MSG2(verb,
00695 tab1 << "integral group " << g << " needs cofacets");
00696 break;
00697 default:
00698 TEUCHOS_TEST_FOR_EXCEPT(1);
00699 }
00700 }
00701
00702 if (allTermsNeedCofacets)
00703 {
00704 cellSpec = AllTermsNeedCofacets;
00705 }
00706 else if (noTermsNeedCofacets)
00707 {
00708 cellSpec = NoTermsNeedCofacets;
00709 }
00710
00711 if (!isMaximalCell && ee->maxDiffOrderOnDiscreteFunctions() > 0)
00712 {
00713 Tabs tab1;
00714 SUNDANCE_MSG2(verb, tab1
00715 << "(*) discrete functions will require cofacet-based transformations");
00716 if (cellSpec==NoTermsNeedCofacets)
00717 {
00718 cellSpec = SomeTermsNeedCofacets;
00719 }
00720 }
00721
00722 SUNDANCE_MSG2(verb, tab << "found: " << cellSpec);
00723
00724 return cellSpec;
00725 }
00726
00727
00728 void Assembler::configureVector(Array<Vector<double> >& b) const
00729 {
00730
00731 TimeMonitor timer(configTimer());
00732
00733 Tabs tab0;
00734 int verb = eqn_->maxWatchFlagSetting("vector config");
00735
00736 SUNDANCE_MSG1(verb, tab0 << "in Assembler::configureVector()");
00737
00738
00739 Array<VectorSpace<double> > vs(eqn_->numVarBlocks());
00740 for (int i=0; i<eqn_->numVarBlocks(); i++)
00741 {
00742 vs[i] = privateRowSpace_[i]->vecSpace();
00743 }
00744 VectorSpace<double> rowSpace;
00745
00746
00747
00748 if (eqn_->numVarBlocks() > 1)
00749 {
00750 rowSpace = Playa::blockSpace(vs);
00751 }
00752 else
00753 {
00754 rowSpace = vs[0];
00755 }
00756
00757
00758 for (int i=0; i<b.size(); i++)
00759 {
00760
00761
00762 b[i] = rowSpace.createMember();
00763
00764
00765 if (!partitionBCs_ && eqn_->numVarBlocks() > 1)
00766 {
00767
00768 Vector<double> vecBlock;
00769 for (int br=0; br<eqn_->numVarBlocks(); br++)
00770 {
00771 configureVectorBlock(br, vecBlock);
00772 b[i].setBlock(br, vecBlock);
00773 }
00774 }
00775 else
00776 {
00777
00778 if (!partitionBCs_)
00779 {
00780 Playa::LoadableVector<double>* lv
00781 = dynamic_cast<Playa::LoadableVector<double>* >(b[i].ptr().get());
00782
00783 TEUCHOS_TEST_FOR_EXCEPTION(lv == 0, std::runtime_error,
00784 "vector is not loadable in Assembler::configureVector()");
00785 }
00786 else
00787 {
00788 }
00789 }
00790 }
00791 numConfiguredColumns_ = b.size();
00792 }
00793
00794 void Assembler::configureVectorBlock(int br, Vector<double>& b) const
00795 {
00796 Tabs tab0;
00797 int verb = eqn_->maxWatchFlagSetting("vector config");
00798 SUNDANCE_MSG2(verb, tab0 << "in Assembler::configureVectorBlock()");
00799 VectorSpace<double> vecSpace = privateRowSpace_[br]->vecSpace();
00800
00801 b = vecSpace.createMember();
00802
00803 if (!partitionBCs_)
00804 {
00805 Playa::LoadableVector<double>* lv
00806 = dynamic_cast<Playa::LoadableVector<double>* >(b.ptr().get());
00807
00808 TEUCHOS_TEST_FOR_EXCEPTION(lv == 0, std::runtime_error,
00809 "vector block is not loadable "
00810 "in Assembler::configureVectorBlock()");
00811 }
00812 }
00813
00814
00815 void Assembler::configureMatrix(LinearOperator<double>& A,
00816 Array<Vector<double> >& b) const
00817 {
00818 TimeMonitor timer(configTimer());
00819 int verb = eqn_->maxWatchFlagSetting("matrix config");
00820 Tabs tab;
00821 SUNDANCE_MSG1(verb, tab << "in Assembler::configureMatrix()");
00822
00823
00824 SUNDANCE_MSG1(verb, tab << "before config: A=" << A.description());
00825 if (matNeedsConfiguration_)
00826 {
00827 Tabs tab0;
00828
00829 int nRowBlocks = rowMap_.size();
00830 int nColBlocks = colMap_.size();
00831 Array<Array<int> > isNonzero = findNonzeroBlocks();
00832
00833 if (nRowBlocks==1 && nColBlocks==1)
00834 {
00835 configureMatrixBlock(0,0,A);
00836 }
00837 else
00838 {
00839 A = makeBlockOperator(solnVecSpace(), rowVecSpace());
00840 for (int br=0; br<nRowBlocks; br++)
00841 {
00842 for (int bc=0; bc<nColBlocks; bc++)
00843 {
00844 if (isNonzero[br][bc])
00845 {
00846 LinearOperator<double> matBlock;
00847 configureMatrixBlock(br, bc, matBlock);
00848 A.setBlock(br, bc, matBlock);
00849 }
00850 }
00851 }
00852 A.endBlockFill();
00853 }
00854 matNeedsConfiguration_ = false;
00855 }
00856 else
00857 {
00858 Tabs tab0;
00859 SUNDANCE_MSG1(verb,
00860 tab0 << "Assembler::configureMatrix() not needed, proceeding to configure vector");
00861 }
00862 SUNDANCE_MSG1(verb, tab << "after config: A=" << A.description());
00863 configureVector(b);
00864 }
00865
00866 void Assembler::configureMatrixBlock(int br, int bc,
00867 LinearOperator<double>& A) const
00868 {
00869 Tabs tab;
00870 TimeMonitor timer(configTimer());
00871 int verb = eqn_->maxWatchFlagSetting("matrix config");
00872
00873 SUNDANCE_MSG1(verb, tab << "in configureMatrixBlock()");
00874
00875 SUNDANCE_MSG2(verb, tab << "num rows = " << rowMap()[br]->numDOFs());
00876
00877 SUNDANCE_MSG2(verb, tab << "num cols = " << colMap()[bc]->numDOFs());
00878
00879 VectorSpace<double> rowSpace = privateRowSpace_[br]->vecSpace();
00880 VectorSpace<double> colSpace = privateColSpace_[bc]->vecSpace();
00881
00882 RCP<MatrixFactory<double> > matFactory ;
00883
00884 matFactory = rowVecType_[br].createMatrixFactory(colSpace, rowSpace);
00885
00886 IncrementallyConfigurableMatrixFactory* icmf
00887 = dynamic_cast<IncrementallyConfigurableMatrixFactory*>(matFactory.get());
00888
00889 CollectivelyConfigurableMatrixFactory* ccmf
00890 = dynamic_cast<CollectivelyConfigurableMatrixFactory*>(matFactory.get());
00891
00892 TEUCHOS_TEST_FOR_EXCEPTION(ccmf==0 && icmf==0, std::runtime_error,
00893 "Neither incremental nor collective matrix structuring "
00894 "appears to be available");
00895
00896
00897
00898
00899 if (false )
00900 {
00901 Tabs tab1;
00902 SUNDANCE_MSG2(verb, tab1 << "Assembler: doing collective matrix structuring...");
00903 Array<int> graphData;
00904 Array<int> nnzPerRow;
00905 Array<int> rowPtrs;
00906
00907 using Teuchos::createVector;
00908
00909 getGraph(br, bc, graphData, rowPtrs, nnzPerRow);
00910 ccmf->configure(lowestRow_[br], createVector(rowPtrs), createVector(nnzPerRow), createVector(graphData));
00911 }
00912 else
00913 {
00914 Tabs tab1;
00915 SUNDANCE_MSG2(verb, tab1 << "Assembler: doing incremental matrix structuring...");
00916 incrementalGetGraph(br, bc, icmf);
00917 {
00918 TimeMonitor timer1(matFinalizeTimer());
00919 icmf->finalize();
00920 }
00921 }
00922
00923 SUNDANCE_MSG3(verb, tab << "Assembler: allocating matrix...");
00924 {
00925 TimeMonitor timer1(matAllocTimer());
00926 A = matFactory->createMatrix();
00927 }
00928 }
00929
00930 Playa::LinearOperator<double> Assembler::allocateMatrix() const
00931 {
00932 LinearOperator<double> A;
00933 Array<Vector<double> > b;
00934 configureMatrix(A, b);
00935 return A;
00936 }
00937
00938
00939
00940
00941
00942
00943 void Assembler::displayEvaluationResults(
00944 const EvalContext& context,
00945 const EvaluatableExpr* evalExpr,
00946 const Array<double>& constantCoeffs,
00947 const Array<RCP<EvalVector> >& vectorCoeffs) const
00948 {
00949 Tabs tab;
00950 FancyOStream& os = Out::os();
00951
00952 os << tab << "evaluation results: " << std::endl;
00953
00954 const RCP<SparsitySuperset>& sparsity
00955 = evalExpr->sparsitySuperset(context);
00956
00957 sparsity->print(os, vectorCoeffs, constantCoeffs);
00958 }
00959
00960
00961
00962 void Assembler::assemblyLoop(const ComputationType& compType,
00963 RCP<AssemblyKernelBase> kernel) const
00964 {
00965 Tabs tab;
00966 int verb = 0;
00967 if (eqn_->hasActiveWatchFlag()) verb = max(eqn_->maxWatchFlagSetting("assembly loop"), 1);
00968
00969
00970 SUNDANCE_BANNER1(verb, tab, "Assembly loop");
00971
00972 SUNDANCE_MSG1(verb, tab << "computation type is " << compType);
00973
00974
00975
00976
00977 SUNDANCE_MSG2(verb, tab << "work set size is " << workSetSize());
00978 RCP<Array<int> > workSet = rcp(new Array<int>());
00979 workSet->reserve(workSetSize());
00980
00981
00982
00983
00984 RCP<Array<int> > isLocalFlag = rcp(new Array<int>());
00985 isLocalFlag->reserve(workSetSize());
00986
00987
00988
00989
00990
00991
00992 Array<RCP<EvalVector> > vectorCoeffs;
00993 Array<double> constantCoeffs;
00994
00995
00996 RCP<Array<double> > localValues = rcp(new Array<double>());
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006 const Array<EvalContext>& contexts = contexts_.get(compType);
01007
01008 const Array<Array<RCP<IntegralGroup> > >& groups = groups_.get(compType);
01009
01010 const Array<const EvaluatableExpr*>& evalExprs
01011 = evalExprs_.get(compType);
01012
01013
01014 const Array<int>& skipRqc = skipRqc_.get(compType);
01015
01016
01017
01018
01019
01020 SUNDANCE_MSG1(verb,
01021 tab << "---------- outer assembly loop over subregions");
01022
01023
01024
01025
01026
01027 Array< Array<RCP<AssemblyTransformationBuilder> > > transformations;
01028
01029
01030 transformations.resize(rqc_.size());
01031 for (int r=0; r<rqc_.size(); r++)
01032 {
01033 transformations[r].resize(groups[r].size());
01034 for (int g=0; g<groups[r].size(); g++)
01035 {
01036 const RCP<IntegralGroup>& group = groups[r][g];
01037
01038
01039 transformations[r][g] = rcp(new AssemblyTransformationBuilder( group , rowMap_ , colMap_ , mesh_));
01040 }
01041 }
01042
01043
01044
01045
01046 int oldKernelVerb = kernel->verb();
01047
01048 TEUCHOS_TEST_FOR_EXCEPT(rqc_.size() != evalExprs.size());
01049
01050
01051 for (int r=0; r<rqc_.size(); r++)
01052 {
01053 Tabs tab0;
01054
01055
01056 int rqcVerb = 0;
01057 int evalVerb = 0;
01058 int evalMedVerb = 0;
01059 int dfEvalVerb = 0;
01060 int fillVerb = 0;
01061
01062
01063 if (rqc_[r].watch().isActive())
01064 {
01065 Tabs tab01;
01066 rqcVerb=verb;
01067 evalVerb = rqc_[r].watch().param("evaluation");
01068 evalMedVerb = rqc_[r].watch().param("eval mediator");
01069 dfEvalVerb = rqc_[r].watch().param("discrete function evaluation");
01070 fillVerb = rqc_[r].watch().param("fill");
01071
01072 SUNDANCE_MSG1(rqcVerb, tab0 << std::endl
01073 << tab0 << "-------------"
01074 << std::endl << tab0 << " doing watched subregion r=" << r << " of "
01075 << rqc_.size() << ", rqc="
01076 << rqc_[r]);
01077 if (skipRqc[r])
01078 {
01079 SUNDANCE_MSG2(rqcVerb, tab01 << "this rqc is not needed in comp type = " << compType);
01080 }
01081 else
01082 {
01083 SUNDANCE_MSG2(rqcVerb, tab01 << "expr is " << evalExprs[r]->toString());
01084 SUNDANCE_MSG2(rqcVerb, tab01 << "isBC= " << isBCRqc_[r]);
01085 }
01086 }
01087 else
01088 {
01089 SUNDANCE_MSG1(verb, tab0 << "unwatched region r=" << r << " of " << rqc_.size());
01090 }
01091 Tabs tab01;
01092
01093 kernel->setVerb(fillVerb);
01094
01095
01096
01097
01098 if ((!isBCRqc_[r] && eqn_->skipRqc(compType, rqc_[r]))
01099 || (isBCRqc_[r] && eqn_->skipBCRqc(compType, rqc_[r])))
01100 {
01101 Tabs tab012;
01102 SUNDANCE_MSG2(rqcVerb, tab012 << "nothing to do for comp type "
01103 << compType);
01104 continue;
01105 }
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116 evalMgr_->setMediator(mediators_[r]);
01117 mediators_[r]->setVerb(evalMedVerb, dfEvalVerb);
01118
01119
01120
01121
01122 evalMgr_->setRegion(contexts_.get(compType)[r]);
01123
01124
01125 CellFilter filter = rqc_[r].domain();
01126
01127
01128
01129
01130
01131 CellSet cells = filter.getCells(mesh_);
01132
01133 int cellDim = filter.dimension(mesh_);
01134
01135
01136
01137
01138
01139
01140 CellType cellType = mesh_.cellType(cellDim);
01141
01142
01143
01144
01145 CellType maxCellType = mesh_.cellType(mesh_.spatialDim());
01146
01147 SUNDANCE_MSG2(rqcVerb, tab01 << "cell type = " << cellType
01148 << ", cellDim = " << cellDim
01149 << ", max cell type = " << maxCellType
01150 << ", max cell dim = " << mesh_.spatialDim());
01151
01152
01153
01154
01155 IntegrationCellSpecifier intCellSpec
01156 = rqcRequiresMaximalCofacets_.get(compType)[r];
01157 SUNDANCE_MSG2(rqcVerb, tab01
01158 << "whether we need to refer to maximal cofacets: " << intCellSpec);
01159
01160
01161
01162 const Array<Set<int> >& requiredVars = eqn_->reducedVarsOnRegion(filter);
01163 const Array<Set<int> >& requiredUnks = eqn_->reducedUnksOnRegion(filter);
01164
01165
01166
01167
01168 SUNDANCE_MSG2(rqcVerb, tab01 << "setting integration specifier in mediator");
01169 mediators_[r]->setIntegrationSpec(intCellSpec);
01170
01171
01172
01173
01174
01175
01176
01177 SUNDANCE_MSG2(rqcVerb, tab01 << "setting cell type=" << cellType << " in mediator");
01178 mediators_[r]->setCellType(cellType, maxCellType, isInternalBdry_[r]);
01179
01180
01181
01182 const Evaluator* evaluator
01183 = evalExprs[r]->evaluator(contexts[r]).get();
01184
01185
01186
01187
01188 CellIterator iter=cells.begin();
01189 int workSetCounter = 0;
01190 int myRank = mesh_.comm().getRank();
01191
01192 SUNDANCE_MSG2(rqcVerb, tab01 << "----- looping over worksets");
01193 while (iter != cells.end())
01194 {
01195 Tabs tab1;
01196
01197
01198
01199
01200
01201
01202 workSet->resize(0);
01203 isLocalFlag->resize(0);
01204 for (int c=0; c<workSetSize() && iter != cells.end(); c++, iter++)
01205 {
01206 workSet->append(*iter);
01207
01208
01209 isLocalFlag->append(myRank==mesh_.ownerProcID(cellDim, *iter));
01210 }
01211
01212 SUNDANCE_MSG2(rqcVerb,
01213 tab1 << "doing work set " << workSetCounter
01214 << " consisting of " << workSet->size() << " cells");
01215 {
01216 Tabs tab2;
01217 SUNDANCE_MSG4(rqcVerb, tab2 << "cells are " << *workSet);
01218 }
01219 workSetCounter++;
01220
01221
01222 evalMgr_->setVerb(evalVerb);
01223
01224
01225
01226
01227 mediators_[r]->setCellBatch(workSet);
01228
01229
01230
01231
01232
01233
01234
01235 const CellJacobianBatch& JVol = mediators_[r]->JVol();
01236 const CellJacobianBatch& JTrans = mediators_[r]->JTrans();
01237
01238
01239
01240
01241 const Array<int>& facetIndices = mediators_[r]->facetIndices();
01242 if (facetIndices.size() > 0)
01243 {
01244 Tabs tab2;
01245 SUNDANCE_MSG2(rqcVerb, tab2 << "facet indices are " << facetIndices);
01246 }
01247
01248
01249
01250
01251 kernel->prepareForWorkSet(
01252 requiredVars,
01253 requiredUnks,
01254 mediators_[r]);
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266 evaluator->resetNumCalls();
01267 SUNDANCE_MSG2(rqcVerb, tab1
01268 << "====== evaluating coefficient expressions") ;
01269 try
01270 {
01271 evalExprs[r]->evaluate(*evalMgr_, constantCoeffs, vectorCoeffs);
01272 }
01273 catch(std::exception& exc)
01274 {
01275 Tabs tabX;
01276 SUNDANCE_BANNER1(10, tabX,
01277 "DETECTED EXCEPTION DURING EXPR EVALUATION CALL IN ASSEMBLY LOOP");
01278 Tabs tabX1;
01279 SUNDANCE_MSG1(10, tabX1 << "While working on RQC="
01280 << rqc_[r]);
01281 SUNDANCE_MSG1(10, tabX1 << "While evaluating expr="
01282 << evalExprs[r]->toString());
01283 throw (std::runtime_error(exc.what()));
01284 }
01285
01286
01287 SUNDANCE_MSG2(rqcVerb, tab1 << "====== done evaluating expressions") ;
01288 if (evalVerb > 3)
01289 {
01290 displayEvaluationResults(contexts[r], evalExprs[r], constantCoeffs,
01291 vectorCoeffs);
01292 }
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306 ElementIntegral::invalidateTransformationMatrices();
01307
01308
01309 SUNDANCE_MSG2(rqcVerb, tab1 << "----- looping over integral groups");
01310 for (int g=0; g<groups[r].size(); g++)
01311 {
01312 Tabs tab2;
01313 SUNDANCE_MSG2(rqcVerb, tab2 << std::endl << tab2
01314 << "--- evaluating integral group g=" << g << " of "
01315 << groups[r].size() );
01316
01317
01318
01319 const RCP<IntegralGroup>& group = groups[r][g];
01320 if (!group->evaluate(JTrans, JVol, *isLocalFlag, facetIndices, workSet,
01321 vectorCoeffs, constantCoeffs, localValues)) continue;
01322
01323
01324
01325
01326
01327 transformations[r][g]->applyTransformsToAssembly(
01328 g , (localValues->size() / workSet->size()),
01329 cellType, cellDim , maxCellType,
01330 JTrans, JVol, facetIndices, workSet, localValues);
01331
01332
01333
01334
01335
01336
01337
01338
01339 {
01340 TimeMonitor fillTM(fillTimer());
01341 kernel->fill(isBCRqc_[r], *group, localValues);
01342 }
01343 }
01344 SUNDANCE_MSG2(rqcVerb, tab1 << "----- done looping over integral groups");
01345 }
01346 SUNDANCE_MSG2(rqcVerb, tab0 << "----- done looping over worksets");
01347
01348 kernel->setVerb(oldKernelVerb);
01349 SUNDANCE_MSG1(verb, tab0 << "----- done rqc");
01350 }
01351 SUNDANCE_MSG1(verb, tab << "----- done looping over rqcs");
01352
01353
01354
01355 SUNDANCE_MSG2(verb, tab << "doing post-loop processing");
01356 kernel->postLoopFinalization();
01357 SUNDANCE_BANNER1(verb, tab, "done assembly loop");
01358
01359
01360 }
01361
01362
01363
01364
01365
01366 void Assembler::assemble(LinearOperator<double>& A,
01367 Array<Vector<double> >& mv) const
01368 {
01369 TimeMonitor timer(assemblyTimer());
01370 Tabs tab;
01371 int verb = 0;
01372 if (eqn_->hasActiveWatchFlag()) verb = max(verb, 1);
01373
01374 SUNDANCE_BANNER1(verb, tab, "Assembling matrix and vector");
01375
01376 TEUCHOS_TEST_FOR_EXCEPTION(!contexts_.containsKey(MatrixAndVector),
01377 std::runtime_error,
01378 "Assembler::assemble(A, b) called for an assembler that "
01379 "does not support matrix/vector assembly");
01380
01381 configureMatrix(A, mv);
01382
01383 RCP<AssemblyKernelBase> kernel
01384 = rcp(new MatrixVectorAssemblyKernel(
01385 rowMap_, isBCRow_, lowestRow_,
01386 colMap_, isBCCol_, lowestCol_,
01387 A, mv, partitionBCs_,
01388 0));
01389
01390 assemblyLoop(MatrixAndVector, kernel);
01391
01392 SUNDANCE_MSG1(verb, tab << "matrix=" << A);
01393 if (verb>0) A.print(Out::os());
01394 SUNDANCE_MSG1(verb, tab << "vectors=" << mv);
01395 for (int i=0; i<mv.size(); i++)
01396 {
01397 SUNDANCE_MSG1(verb, tab << "vectors #" << i);
01398 if (verb>0) mv[i].print(Out::os());
01399 }
01400
01401 SUNDANCE_MSG1(verb, tab << "Assembler: done assembling matrix and vector");
01402 }
01403
01404
01405
01406 void Assembler::assembleSensitivities(LinearOperator<double>& A,
01407 Array<Vector<double> >& mv) const
01408 {
01409 TimeMonitor timer(assemblyTimer());
01410 Tabs tab;
01411 int verb = 0;
01412 if (eqn_->hasActiveWatchFlag()) verb = max(verb, 1);
01413
01414 SUNDANCE_BANNER1(verb, tab, "Assembling matrix and sensitivity vector");
01415
01416 TEUCHOS_TEST_FOR_EXCEPTION(!contexts_.containsKey(Sensitivities),
01417 std::runtime_error,
01418 "Assembler::assembleSensitivities(A, b) called for an assembler that "
01419 "does not support sensitivity assembly");
01420
01421 configureMatrix(A, mv);
01422
01423
01424 RCP<AssemblyKernelBase> kernel
01425 = rcp(new MatrixVectorAssemblyKernel(
01426 rowMap_, isBCRow_, lowestRow_,
01427 colMap_, isBCCol_, lowestCol_,
01428 A, mv, partitionBCs_,
01429 0));
01430
01431 assemblyLoop(Sensitivities, kernel);
01432 SUNDANCE_MSG1(verb, tab << "Assembler: done assembling matrix and sensitivity vector");
01433 }
01434
01435
01436
01437
01438 void Assembler::assemble(Array<Vector<double> >& mv) const
01439 {
01440
01441 Tabs tab;
01442
01443 TimeMonitor timer(assemblyTimer());
01444
01445
01446 int verb = 0;
01447 if (eqn_->hasActiveWatchFlag()) verb = max(verb, 1);
01448
01449 SUNDANCE_BANNER1(verb, tab, "Assembling vector");
01450
01451
01452 TEUCHOS_TEST_FOR_EXCEPTION(!contexts_.containsKey(VectorOnly),
01453 std::runtime_error,
01454 "Assembler::assemble(b) called for an assembler that "
01455 "does not support vector-only assembly");
01456
01457
01458 configureVector(mv);
01459
01460
01461 RCP<AssemblyKernelBase> kernel
01462 = rcp(new VectorAssemblyKernel(
01463 rowMap_, isBCRow_, lowestRow_,
01464 mv, partitionBCs_, 0));
01465
01466 assemblyLoop(VectorOnly, kernel);
01467
01468 SUNDANCE_MSG1(verb, tab << "Assembler: done assembling vector");
01469 }
01470
01471
01472
01473
01474 void Assembler::evaluate(double& value, Array<Vector<double> >& gradient) const
01475 {
01476 Tabs tab;
01477 TimeMonitor timer(assemblyTimer());
01478 int verb = 0;
01479 if (eqn_->hasActiveWatchFlag()) verb = max(verb, 1);
01480
01481 SUNDANCE_BANNER1(verb, tab, "Computing functional and gradient");
01482
01483 TEUCHOS_TEST_FOR_EXCEPTION(!contexts_.containsKey(FunctionalAndGradient),
01484 std::runtime_error,
01485 "Assembler::evaluate(f,df) called for an assembler that "
01486 "does not support value/gradient assembly");
01487
01488 configureVector(gradient);
01489
01490 value = 0.0;
01491
01492 RCP<AssemblyKernelBase> kernel
01493 = rcp(new FunctionalGradientAssemblyKernel(
01494 mesh_.comm(),
01495 rowMap_, isBCRow_, lowestRow_,
01496 gradient, partitionBCs_, &value, verb));
01497
01498 assemblyLoop(FunctionalAndGradient, kernel);
01499
01500 SUNDANCE_BANNER1(verb, tab, "Done computing functional and gradient");
01501
01502 }
01503
01504
01505
01506
01507
01508
01509 void Assembler::evaluate(double& value) const
01510 {
01511 Tabs tab;
01512 TimeMonitor timer(assemblyTimer());
01513 int verb = 0;
01514 if (eqn_->hasActiveWatchFlag()) verb = max(verb, 1);
01515
01516 SUNDANCE_BANNER1(verb, tab, "Computing functional");
01517
01518 TEUCHOS_TEST_FOR_EXCEPTION(!contexts_.containsKey(FunctionalOnly),
01519 std::runtime_error,
01520 "Assembler::evaluate(f) called for an assembler that "
01521 "does not support functional evaluation");
01522
01523 value = 0.0;
01524
01525 RCP<AssemblyKernelBase> kernel
01526 = rcp(new FunctionalAssemblyKernel(mesh_.comm(), &value, verb));
01527
01528 assemblyLoop(FunctionalOnly, kernel);
01529
01530 SUNDANCE_BANNER1(verb, tab, "Done computing functional");
01531 }
01532
01533
01534
01535
01536
01537
01538
01539 void Assembler::getGraph(int br, int bc,
01540 Array<int>& graphData,
01541 Array<int>& rowPtrs,
01542 Array<int>& nnzPerRow) const
01543 {
01544 TimeMonitor timer(graphBuildTimer());
01545 Tabs tab;
01546 int verb = eqn_->maxWatchFlagSetting("matrix config");
01547
01548
01549 RCP<Array<int> > workSet = rcp(new Array<int>());
01550 workSet->reserve(workSetSize());
01551
01552 RCP<Array<Array<int> > > testLocalDOFs
01553 = rcp(new Array<Array<int> >());
01554
01555 RCP<Array<Array<int> > > unkLocalDOFs
01556 = rcp(new Array<Array<int> >());
01557
01558 SUNDANCE_MSG3(verb, tab << "Creating graph: there are "
01559 << rowMap_[br]->numLocalDOFs()
01560 << " local equations");
01561
01562
01563 Array<Set<int> > tmpGraph;
01564 tmpGraph.resize(rowMap_[br]->numLocalDOFs());
01565
01566 {
01567 TimeMonitor timer2(colSearchTimer());
01568 for (int d=0; d<eqn_->numRegions(); d++)
01569 {
01570 Tabs tab0;
01571 CellFilter domain = eqn_->region(d);
01572 const Array<Set<int> >& requiredVars = eqn_->reducedVarsOnRegion(domain);
01573 const Array<Set<int> >& requiredUnks = eqn_->reducedUnksOnRegion(domain);
01574 SUNDANCE_MSG2(verb,
01575 tab0 << "cell set " << domain
01576 << " isBCRegion=" << eqn_->isBCRegion(d));
01577 int dim = domain.dimension(mesh_);
01578 CellSet cells = domain.getCells(mesh_);
01579
01580 RCP<Set<OrderedPair<int, int> > > pairs ;
01581 if (eqn_->hasVarUnkPairs(domain)) pairs = eqn_->varUnkPairs(domain);
01582
01583 SUNDANCE_OUT(verb > 2 && pairs.get() != 0,
01584 tab0 << "non-BC pairs = "
01585 << *pairs);
01586
01587 RCP<Set<OrderedPair<int, int> > > bcPairs ;
01588 if (eqn_->isBCRegion(d))
01589 {
01590 if (eqn_->hasBCVarUnkPairs(domain))
01591 {
01592 bcPairs = eqn_->bcVarUnkPairs(domain);
01593 SUNDANCE_OUT(verb > 2, tab0 << "BC pairs = "
01594 << *bcPairs);
01595 }
01596 }
01597 Array<Set<int> > unksForTestsSet(eqn_->numVarIDs(bc));
01598 Array<Set<int> > bcUnksForTestsSet(eqn_->numVarIDs(bc));
01599
01600 Set<OrderedPair<int, int> >::const_iterator i;
01601
01602 if (pairs.get() != 0)
01603 {
01604 for (i=pairs->begin(); i!=pairs->end(); i++)
01605 {
01606 const OrderedPair<int, int>& p = *i;
01607 int t = p.first();
01608 int u = p.second();
01609
01610 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasVarID(t), std::logic_error,
01611 "Test function ID " << t << " does not appear "
01612 "in equation set");
01613 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasUnkID(u), std::logic_error,
01614 "Unk function ID " << u << " does not appear "
01615 "in equation set");
01616
01617 if (eqn_->blockForVarID(t) != br) continue;
01618 if (eqn_->blockForUnkID(u) != bc) continue;
01619
01620 unksForTestsSet[eqn_->reducedVarID(t)].put(eqn_->reducedUnkID(u));
01621 }
01622 }
01623 if (bcPairs.get() != 0)
01624 {
01625 for (i=bcPairs->begin(); i!=bcPairs->end(); i++)
01626 {
01627 const OrderedPair<int, int>& p = *i;
01628 int t = p.first();
01629 int u = p.second();
01630
01631 if (eqn_->blockForVarID(t) != br) continue;
01632 if (eqn_->blockForUnkID(u) != bc) continue;
01633
01634 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasVarID(t), std::logic_error,
01635 "Test function ID " << t << " does not appear "
01636 "in equation set");
01637 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasUnkID(u), std::logic_error,
01638 "Unk function ID " << u << " does not appear "
01639 "in equation set");
01640 bcUnksForTestsSet[eqn_->reducedVarID(t)].put(eqn_->reducedUnkID(u));
01641 }
01642 }
01643
01644 Array<Array<int> > unksForTests(unksForTestsSet.size());
01645 Array<Array<int> > bcUnksForTests(bcUnksForTestsSet.size());
01646
01647 for (int t=0; t<unksForTests.size(); t++)
01648 {
01649 unksForTests[t] = unksForTestsSet[t].elements();
01650 bcUnksForTests[t] = bcUnksForTestsSet[t].elements();
01651 }
01652
01653 Array<int> numTestNodes;
01654 Array<int> numUnkNodes;
01655
01656
01657
01658 int highestRow = lowestRow_[br] + rowMap_[br]->numLocalDOFs();
01659
01660 int nt = eqn_->numVarIDs(br);
01661
01662 CellIterator iter=cells.begin();
01663 while (iter != cells.end())
01664 {
01665
01666 workSet->resize(0);
01667 for (int c=0; c<workSetSize() && iter != cells.end(); c++, iter++)
01668 {
01669 workSet->append(*iter);
01670 }
01671
01672 int nCells = workSet->size();
01673
01674 RCP<const MapStructure> colMapStruct;
01675 RCP<const MapStructure> rowMapStruct
01676 = rowMap_[br]->getDOFsForCellBatch(dim, *workSet,
01677 requiredVars[br], *testLocalDOFs,
01678 numTestNodes, verb);
01679 if (rowMap_[br].get()==colMap_[bc].get())
01680 {
01681 unkLocalDOFs = testLocalDOFs;
01682 numUnkNodes = numTestNodes;
01683 colMapStruct = rowMapStruct;
01684 }
01685 else
01686 {
01687 colMapStruct = colMap_[br]->getDOFsForCellBatch(dim, *workSet,
01688 requiredUnks[bc],
01689 *unkLocalDOFs, numUnkNodes, verb);
01690 }
01691
01692 if (pairs.get() != 0)
01693 {
01694 for (int c=0; c<nCells; c++)
01695 {
01696 for (int t=0; t<nt; t++)
01697 {
01698 int tChunk = rowMapStruct->chunkForFuncID(t);
01699 int nTestFuncs = rowMapStruct->numFuncs(tChunk);
01700 int testFuncIndex = rowMapStruct->indexForFuncID(t);
01701 int nTestNodes = numTestNodes[tChunk];
01702 const Array<int>& testDOFs = (*testLocalDOFs)[tChunk];
01703 for (int uit=0; uit<unksForTests[t].size(); uit++)
01704 {
01705 Tabs tab2;
01706 int u = unksForTests[t][uit];
01707 int uChunk = colMapStruct->chunkForFuncID(u);
01708 int nUnkFuncs = colMapStruct->numFuncs(uChunk);
01709 int unkFuncIndex = colMapStruct->indexForFuncID(u);
01710 const Array<int>& unkDOFs = (*unkLocalDOFs)[uChunk];
01711 int nUnkNodes = numUnkNodes[uChunk];
01712 for (int n=0; n<nTestNodes; n++)
01713 {
01714 int row
01715 = testDOFs[(c*nTestFuncs + testFuncIndex)*nTestNodes + n];
01716 if (row < lowestRow_[br] || row >= highestRow
01717 || (*(isBCRow_[br]))[row-lowestRow_[br]]) continue;
01718 Set<int>& colSet = tmpGraph[row-lowestRow_[br]];
01719 for (int m=0; m<nUnkNodes; m++)
01720 {
01721 int col
01722 = unkDOFs[(c*nUnkFuncs + unkFuncIndex)*nUnkNodes + m];
01723 colSet.put(col);
01724 }
01725 }
01726 }
01727 }
01728 }
01729 }
01730 if (bcPairs.get() != 0)
01731 {
01732 for (int c=0; c<nCells; c++)
01733 {
01734 for (int t=0; t<nt; t++)
01735 {
01736 int tChunk = rowMapStruct->chunkForFuncID(t);
01737 int nTestFuncs = rowMapStruct->numFuncs(tChunk);
01738 int testFuncIndex = rowMapStruct->indexForFuncID(t);
01739 int nTestNodes = numTestNodes[tChunk];
01740 const Array<int>& testDOFs = (*testLocalDOFs)[tChunk];
01741 for (int uit=0; uit<bcUnksForTests[t].size(); uit++)
01742 {
01743 Tabs tab2;
01744 int u = bcUnksForTests[t][uit];
01745 int uChunk = colMapStruct->chunkForFuncID(u);
01746 int nUnkFuncs = colMapStruct->numFuncs(uChunk);
01747 int unkFuncIndex = colMapStruct->indexForFuncID(u);
01748 const Array<int>& unkDOFs = (*unkLocalDOFs)[uChunk];
01749 int nUnkNodes = numUnkNodes[uChunk];
01750 for (int n=0; n<nTestNodes; n++)
01751 {
01752 int row
01753 = testDOFs[(c*nTestFuncs + testFuncIndex)*nTestNodes + n];
01754 if (row < lowestRow_[br] || row >= highestRow
01755 || !(*(isBCRow_[br]))[row-lowestRow_[br]]) continue;
01756 Set<int>& colSet = tmpGraph[row-lowestRow_[br]];
01757 for (int m=0; m<nUnkNodes; m++)
01758 {
01759 int col
01760 = unkDOFs[(c*nUnkFuncs + unkFuncIndex)*nUnkNodes + m];
01761 colSet.put(col);
01762 }
01763 }
01764 }
01765 }
01766 }
01767 }
01768 }
01769 }
01770 }
01771
01772
01773 {
01774 TimeMonitor t2(graphFlatteningTimer());
01775 int nLocalRows = rowMap_[br]->numLocalDOFs();
01776
01777 int nnz = 0;
01778 rowPtrs.resize(nLocalRows);
01779 nnzPerRow.resize(rowMap_[br]->numLocalDOFs());
01780 for (int i=0; i<nLocalRows; i++)
01781 {
01782 rowPtrs[i] = nnz;
01783 nnzPerRow[i] = tmpGraph[i].size();
01784 nnz += nnzPerRow[i];
01785 }
01786
01787 graphData.resize(nnz);
01788 int* base = &(graphData[0]);
01789 for (int i=0; i<nLocalRows; i++)
01790 {
01791
01792 int* rowBase = base + rowPtrs[i];
01793 const Set<int>& rowSet = tmpGraph[i];
01794 int k = 0;
01795 for (Set<int>::const_iterator
01796 j=rowSet.begin(); j != rowSet.end(); j++, k++)
01797 {
01798 rowBase[k] = *j;
01799 }
01800 }
01801 }
01802
01803 }
01804
01805
01806
01807 void Assembler
01808 ::incrementalGetGraph(int br, int bc,
01809 IncrementallyConfigurableMatrixFactory* icmf) const
01810 {
01811 TimeMonitor timer(graphBuildTimer());
01812 Tabs tab;
01813 int verb = eqn_->maxWatchFlagSetting("matrix config");
01814
01815 RCP<Array<int> > workSet = rcp(new Array<int>());
01816 workSet->reserve(workSetSize());
01817
01818 RCP<Array<Array<int> > > testLocalDOFs
01819 = rcp(new Array<Array<int> >());
01820
01821 RCP<Array<Array<int> > > unkLocalDOFs
01822 = rcp(new Array<Array<int> >());
01823
01824 SUNDANCE_MSG2(verb, tab << "Creating graph: there are "
01825 << rowMap_[br]->numLocalDOFs()
01826 << " local equations");
01827
01828
01829 for (int d=0; d<eqn_->numRegions(); d++)
01830 {
01831 Tabs tab0;
01832 CellFilter domain = eqn_->region(d);
01833 const Array<Set<int> >& requiredVars = eqn_->reducedVarsOnRegion(domain);
01834 const Array<Set<int> >& requiredUnks = eqn_->reducedUnksOnRegion(domain);
01835 Array<int> localVars = requiredVars[br].elements();
01836 Array<int> localUnks = requiredUnks[bc].elements();
01837 SUNDANCE_MSG3(verb,
01838 tab0 << "cell set " << domain
01839 << " isBCRegion=" << eqn_->isBCRegion(d));
01840 int dim = domain.dimension(mesh_);
01841 CellSet cells = domain.getCells(mesh_);
01842
01843 RCP<Set<OrderedPair<int, int> > > pairs ;
01844 if (eqn_->hasVarUnkPairs(domain)) pairs = eqn_->varUnkPairs(domain);
01845
01846 SUNDANCE_OUT(verb > 2 && pairs.get() != 0,
01847 tab0 << "non-BC pairs = "
01848 << *pairs);
01849
01850 RCP<Set<OrderedPair<int, int> > > bcPairs ;
01851 if (eqn_->isBCRegion(d))
01852 {
01853 if (eqn_->hasBCVarUnkPairs(domain))
01854 {
01855 bcPairs = eqn_->bcVarUnkPairs(domain);
01856 SUNDANCE_MSG3(verb, tab0 << "BC pairs = "
01857 << *bcPairs);
01858 }
01859 }
01860 Array<Set<int> > unksForTestsSet(eqn_->numVarIDs(br));
01861 Array<Set<int> > bcUnksForTestsSet(eqn_->numVarIDs(br));
01862
01863 Set<OrderedPair<int, int> >::const_iterator i;
01864
01865 if (pairs.get() != 0)
01866 {
01867 for (i=pairs->begin(); i!=pairs->end(); i++)
01868 {
01869 const OrderedPair<int, int>& p = *i;
01870 int t = p.first();
01871 int u = p.second();
01872
01873 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasVarID(t), std::logic_error,
01874 "Test function ID " << t << " does not appear "
01875 "in equation set");
01876 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasUnkID(u), std::logic_error,
01877 "Unk function ID " << u << " does not appear "
01878 "in equation set");
01879
01880
01881 if (eqn_->blockForVarID(t) != br) continue;
01882 if (eqn_->blockForUnkID(u) != bc) continue;
01883
01884 unksForTestsSet[eqn_->reducedVarID(t)].put(eqn_->reducedUnkID(u));
01885 }
01886 }
01887 if (bcPairs.get() != 0)
01888 {
01889 for (i=bcPairs->begin(); i!=bcPairs->end(); i++)
01890 {
01891 const OrderedPair<int, int>& p = *i;
01892 int t = p.first();
01893 int u = p.second();
01894 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasVarID(t), std::logic_error,
01895 "Test function ID " << t << " does not appear "
01896 "in equation set");
01897 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasUnkID(u), std::logic_error,
01898 "Unk function ID " << u << " does not appear "
01899 "in equation set");
01900
01901 if (eqn_->blockForVarID(t) != br) continue;
01902 if (eqn_->blockForUnkID(u) != bc) continue;
01903
01904 bcUnksForTestsSet[eqn_->reducedVarID(t)].put(eqn_->reducedUnkID(u));
01905 }
01906 }
01907
01908 Array<Array<int> > unksForTests(unksForTestsSet.size());
01909 Array<Array<int> > bcUnksForTests(bcUnksForTestsSet.size());
01910
01911 for (int t=0; t<unksForTests.size(); t++)
01912 {
01913 unksForTests[t] = unksForTestsSet[t].elements();
01914 bcUnksForTests[t] = bcUnksForTestsSet[t].elements();
01915 }
01916
01917 Array<int> numTestNodes;
01918 Array<int> numUnkNodes;
01919
01920 int highestRow = lowestRow_[br] + rowMap_[br]->numLocalDOFs();
01921
01922 CellIterator iter=cells.begin();
01923
01924 while (iter != cells.end())
01925 {
01926
01927 workSet->resize(0);
01928 for (int c=0; c<workSetSize() && iter != cells.end(); c++, iter++)
01929 {
01930 workSet->append(*iter);
01931 }
01932
01933 int nCells = workSet->size();
01934
01935 RCP<const MapStructure> colMapStruct;
01936 RCP<const MapStructure> rowMapStruct
01937 = rowMap_[br]->getDOFsForCellBatch(dim, *workSet,
01938 requiredVars[br],
01939 *testLocalDOFs,
01940 numTestNodes, verb);
01941
01942 if (rowMap_[br].get()==colMap_[bc].get())
01943 {
01944 unkLocalDOFs = testLocalDOFs;
01945 numUnkNodes = numTestNodes;
01946 colMapStruct = rowMapStruct;
01947 }
01948 else
01949 {
01950 colMapStruct = colMap_[bc]->getDOFsForCellBatch(dim, *workSet,
01951 requiredUnks[bc],
01952 *unkLocalDOFs, numUnkNodes, verb);
01953 }
01954
01955
01956 if (pairs.get() != 0)
01957 {
01958 for (int c=0; c<nCells; c++)
01959 {
01960 for (int tIndex=0; tIndex<localVars.size(); tIndex++)
01961 {
01962 int t = localVars[tIndex];
01963 int tChunk = rowMapStruct->chunkForFuncID(t);
01964 int nTestFuncs = rowMapStruct->numFuncs(tChunk);
01965 int testFuncIndex = rowMapStruct->indexForFuncID(t);
01966 int nTestNodes = numTestNodes[tChunk];
01967 const Array<int>& testDOFs = (*testLocalDOFs)[tChunk];
01968 for (int uit=0; uit<unksForTests[t].size(); uit++)
01969 {
01970 Tabs tab2;
01971 int u = unksForTests[t][uit];
01972 int uChunk = colMapStruct->chunkForFuncID(u);
01973 int nUnkFuncs = colMapStruct->numFuncs(uChunk);
01974 int unkFuncIndex = colMapStruct->indexForFuncID(u);
01975 const Array<int>& unkDOFs = (*unkLocalDOFs)[uChunk];
01976 int nUnkNodes = numUnkNodes[uChunk];
01977 for (int n=0; n<nTestNodes; n++)
01978 {
01979 int row
01980 = testDOFs[(c*nTestFuncs + testFuncIndex)*nTestNodes + n];
01981 if (row < lowestRow_[br] || row >= highestRow
01982 || (*(isBCRow_[br]))[row-lowestRow_[br]]) continue;
01983 const int* colPtr = &(unkDOFs[(c*nUnkFuncs + unkFuncIndex)*nUnkNodes]);
01984 icmf->initializeNonzerosInRow(row, nUnkNodes, colPtr);
01985 }
01986 }
01987 }
01988 }
01989 }
01990 if (bcPairs.get() != 0)
01991 {
01992 for (int c=0; c<nCells; c++)
01993 {
01994 for (int tIndex=0; tIndex<localVars.size(); tIndex++)
01995 {
01996 int t = localVars[tIndex];
01997 int tChunk = rowMapStruct->chunkForFuncID(t);
01998 int nTestFuncs = rowMapStruct->numFuncs(tChunk);
01999 int testFuncIndex = rowMapStruct->indexForFuncID(t);
02000 int nTestNodes = numTestNodes[tChunk];
02001 const Array<int>& testDOFs = (*testLocalDOFs)[tChunk];
02002 for (int uit=0; uit<bcUnksForTests[t].size(); uit++)
02003 {
02004 Tabs tab2;
02005 int u = bcUnksForTests[t][uit];
02006 int uChunk = colMapStruct->chunkForFuncID(u);
02007 int nUnkFuncs = colMapStruct->numFuncs(uChunk);
02008 int unkFuncIndex = colMapStruct->indexForFuncID(u);
02009 const Array<int>& unkDOFs = (*unkLocalDOFs)[uChunk];
02010 int nUnkNodes = numUnkNodes[uChunk];
02011 for (int n=0; n<nTestNodes; n++)
02012 {
02013 int row
02014 = testDOFs[(c*nTestFuncs + testFuncIndex)*nTestNodes + n];
02015 if (row < lowestRow_[br] || row >= highestRow
02016 || !(*(isBCRow_[br]))[row-lowestRow_[br]]) continue;
02017
02018 const int* colPtr = &(unkDOFs[(c*nUnkFuncs + unkFuncIndex)*nUnkNodes]);
02019 icmf->initializeNonzerosInRow(row, nUnkNodes, colPtr);
02020 }
02021 }
02022 }
02023 }
02024 }
02025 }
02026 }
02027 }
02028
02029
02030 Array<Array<int> > Assembler::findNonzeroBlocks() const
02031 {
02032 int verb = eqn_->maxWatchFlagSetting("assembler setup");
02033
02034 Array<Array<int> > rtn(eqn_->numVarBlocks());
02035 for (int br=0; br<rtn.size(); br++)
02036 {
02037 rtn[br].resize(eqn_->numUnkBlocks());
02038 for (int bc=0; bc<rtn[br].size(); bc++)
02039 {
02040 rtn[br][bc] = 0 ;
02041 }
02042 }
02043
02044 for (int d=0; d<eqn_->numRegions(); d++)
02045 {
02046 Tabs tab0;
02047 CellFilter domain = eqn_->region(d);
02048 SUNDANCE_MSG3(verb,
02049 tab0 << "cell set " << domain
02050 << " isBCRegion=" << eqn_->isBCRegion(d));
02051
02052 RCP<Set<OrderedPair<int, int> > > pairs ;
02053 if (eqn_->hasVarUnkPairs(domain)) pairs = eqn_->varUnkPairs(domain);
02054
02055 SUNDANCE_OUT(verb > 2 && pairs.get() != 0,
02056 tab0 << "non-BC pairs = "
02057 << *pairs);
02058
02059 RCP<Set<OrderedPair<int, int> > > bcPairs ;
02060 if (eqn_->isBCRegion(d))
02061 {
02062 if (eqn_->hasBCVarUnkPairs(domain))
02063 {
02064 bcPairs = eqn_->bcVarUnkPairs(domain);
02065 SUNDANCE_MSG3(verb, tab0 << "BC pairs = "
02066 << *bcPairs);
02067 }
02068 }
02069
02070 Set<OrderedPair<int, int> >::const_iterator i;
02071
02072 if (pairs.get() != 0)
02073 {
02074 for (i=pairs->begin(); i!=pairs->end(); i++)
02075 {
02076 const OrderedPair<int, int>& p = *i;
02077 int t = p.first();
02078 int u = p.second();
02079
02080 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasVarID(t), std::logic_error,
02081 "Test function ID " << t << " does not appear "
02082 "in equation set");
02083 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasUnkID(u), std::logic_error,
02084 "Unk function ID " << u << " does not appear "
02085 "in equation set");
02086
02087
02088 int br = eqn_->blockForVarID(t);
02089 int bc = eqn_->blockForUnkID(u);
02090 rtn[br][bc] = 1;
02091 }
02092 }
02093 if (bcPairs.get() != 0)
02094 {
02095 for (i=bcPairs->begin(); i!=bcPairs->end(); i++)
02096 {
02097 const OrderedPair<int, int>& p = *i;
02098 int t = p.first();
02099 int u = p.second();
02100 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasVarID(t), std::logic_error,
02101 "Test function ID " << t << " does not appear "
02102 "in equation set");
02103 TEUCHOS_TEST_FOR_EXCEPTION(!eqn_->hasUnkID(u), std::logic_error,
02104 "Unk function ID " << u << " does not appear "
02105 "in equation set");
02106 int br = eqn_->blockForVarID(t);
02107 int bc = eqn_->blockForUnkID(u);
02108 rtn[br][bc] = 1;
02109 }
02110 }
02111 }
02112
02113 return rtn;
02114 }
02115
02116 VectorSpace<double> Assembler::solnVecSpace() const
02117 {
02118 Array<VectorSpace<double> > rtn(eqn_->numUnkBlocks());
02119
02120 for (int i=0; i<rtn.size(); i++)
02121 {
02122 rtn[i] = solutionSpace()[i]->vecSpace();
02123 }
02124
02125 if ((int) rtn.size() == 1)
02126 {
02127 return rtn[0];
02128 }
02129 return blockSpace(rtn);
02130 }
02131
02132
02133 VectorSpace<double> Assembler::rowVecSpace() const
02134 {
02135 Array<VectorSpace<double> > rtn(eqn_->numVarBlocks());
02136
02137 for (int i=0; i<rtn.size(); i++)
02138 {
02139 rtn[i] = rowSpace()[i]->vecSpace();
02140 }
02141
02142 if ((int) rtn.size() == 1)
02143 {
02144 return rtn[0];
02145 }
02146 return blockSpace(rtn);
02147 }
02148
02149
02150
02151 int& Assembler::workSetSize()
02152 {
02153 static int rtn = defaultWorkSetSize();
02154 return rtn;
02155 }
02156
02157 int Assembler::maxWatchFlagSetting(const std::string& name) const
02158 {
02159 return eqnSet()->maxWatchFlagSetting(name);
02160 }
02161
02162
02163
02164 void Assembler::flushConfiguration() const
02165 {
02166 numConfiguredColumns_ = 0;
02167 matNeedsConfiguration_ = true;
02168 mesh_.flushSpecialWeights();
02169 mesh_.flushCurvePoints();
02170
02171
02172 for (int r=0; r<rqc_.size(); r++)
02173 {
02174 const CellFilter filter = rqc_[r].domain();
02175 filter.cfbPtr()->flushCache();
02176 }
02177 }