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 "SundanceFunctionSupportResolver.hpp"
00032 #include "PlayaTabs.hpp"
00033 #include "SundanceSumOfBCs.hpp"
00034 #include "SundanceRegionQuadCombo.hpp"
00035 #include "SundanceOut.hpp"
00036 #include "SundanceUnknownFuncElement.hpp"
00037 #include "SundanceSpectralExpr.hpp"
00038 #include "SundanceUnknownParameterElement.hpp"
00039 #include "SundanceTestFuncElement.hpp"
00040 #include "PlayaExceptions.hpp"
00041
00042 using namespace Sundance;
00043 using namespace Sundance;
00044 using namespace Teuchos;
00045 using namespace std;
00046
00047 FunctionSupportResolver::FunctionSupportResolver(
00048 const Expr& eqns,
00049 const Expr& bcs,
00050 const Array<Expr>& vars,
00051 const Array<Expr>& unks,
00052 const Expr& unkParams,
00053 const Expr& fixedParams,
00054 const Array<Expr>& fixedFields,
00055 bool isVariational)
00056 :
00057 eqns_(eqns),
00058 bcs_(bcs),
00059 integralSum_(dynamic_cast<const SumOfIntegrals*>(eqns.ptr().get())),
00060 bcSum_(0),
00061 varFuncSet_(),
00062 unkFuncSet_(),
00063 unkParamSet_(),
00064 fixedParamSet_(),
00065 regions_(),
00066 regionToIndexMap_(),
00067 varsOnRegions_(),
00068 unksOnRegions_(),
00069 bcVarsOnRegions_(),
00070 bcUnksOnRegions_(),
00071 testToRegionsMap_(),
00072 unkToRegionsMap_(),
00073 varFuncData_(),
00074 unkFuncData_(),
00075 varFuncs_(flattenSpectral(vars)),
00076 unkFuncs_(flattenSpectral(unks)),
00077 fixedFields_(flattenSpectral(fixedFields)),
00078 unkParams_(unkParams),
00079 fixedParams_(fixedParams),
00080 varIDToReducedIDMap_(),
00081 unkIDToReducedIDMap_(),
00082 unkParamIDToReducedUnkParamIDMap_(),
00083 fixedParamIDToReducedFixedParamIDMap_(),
00084 varIDToBlockMap_(),
00085 unkIDToBlockMap_(),
00086 unreducedVarID_(),
00087 unreducedUnkID_(),
00088 unreducedUnkParamID_(),
00089 unreducedFixedParamID_(),
00090 isVariationalProblem_(isVariational)
00091 {
00092 Tabs tab0(0);
00093 Tabs tab1;
00094
00095
00096
00097
00098 TEUCHOS_TEST_FOR_EXCEPTION(eqns.ptr().get()==0, std::runtime_error,
00099 "FunctionSupportResolver ctor detected empty equation set input");
00100
00101 TEUCHOS_TEST_FOR_EXCEPTION(integralSum_==0, std::runtime_error,
00102 "FunctionSupportResolver ctor detected an input equation set that "
00103 "is not in integral form");
00104
00105 bool hasBCs = false;
00106 if (bcs.ptr().get() != 0)
00107 {
00108 bcSum_ = dynamic_cast<const SumOfBCs*>(bcs.ptr().get());
00109 TEUCHOS_TEST_FOR_EXCEPTION(bcSum_==0, std::runtime_error,
00110 "FunctionSupport ctor detected an input Essential "
00111 "BC that is not an EssentialBC object. "
00112 "Object is " << bcs);
00113 hasBCs = true;
00114 }
00115
00116
00117 int verb = 0;
00118 if (integralSum_->hasWatchedTerm() || (hasBCs && bcSum_->hasWatchedTerm()))
00119 {
00120 int v1 = integralSum_->eqnSetSetupVerb();
00121 int v2 = 0;
00122 if (hasBCs) v2 = bcSum_->eqnSetSetupVerb();
00123 verb = std::max(verb, max(v1, v2));
00124 }
00125 SUNDANCE_BANNER1(verb, tab0, "FunctionSupportResolver setup");
00126
00127 if (hasBCs)
00128 {
00129 SUNDANCE_MSG1(verb, tab1 << "Problem has EssentialBCs");
00130 }
00131 else
00132 {
00133 SUNDANCE_MSG1(verb, tab1 << "Problem has no EssentialBCs");
00134 }
00135
00136 SUNDANCE_MSG1(verb, tab1 << "verbosity is " << verb);
00137
00138
00139
00140
00141
00142
00143
00144 bool varsAreTestFunctions = false;
00145 for (int b=0; b<vars.size(); b++)
00146 {
00147 for (int i=0; i<vars[b].size(); i++)
00148 {
00149 const TestFuncElement* t
00150 = dynamic_cast<const TestFuncElement*>(vars[b][i].ptr().get());
00151 TEUCHOS_TEST_FOR_EXCEPTION(t == 0 && varsAreTestFunctions==true,
00152 std::runtime_error,
00153 "variational function list " << vars
00154 << " contains a mix of test and "
00155 "non-test functions");
00156 TEUCHOS_TEST_FOR_EXCEPTION(t != 0 && i>0 && varsAreTestFunctions==false,
00157 std::runtime_error,
00158 "variational function list " << vars
00159 << " contains a mix of test and "
00160 "non-test functions");
00161 if (t!=0) varsAreTestFunctions=true;
00162 }
00163 }
00164
00165 TEUCHOS_TEST_FOR_EXCEPTION(varsAreTestFunctions == true
00166 && isVariationalProblem_,
00167 std::runtime_error,
00168 "test functions given to a variational equation set");
00169
00170 TEUCHOS_TEST_FOR_EXCEPTION(varsAreTestFunctions == false
00171 && !isVariationalProblem_,
00172 std::runtime_error,
00173 "variational functions are unknowns, but equation "
00174 "set is not variational");
00175
00176 if (isVariationalProblem_)
00177 {
00178 SUNDANCE_MSG1(verb, tab1 << "is variational problem");
00179 }
00180 else
00181 {
00182 SUNDANCE_MSG1(verb, tab1 << "is not in variational form");
00183 }
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193 varFuncData_.resize(vars.size());
00194 unkFuncData_.resize(unks.size());
00195 varIDToReducedIDMap_.resize(vars.size());
00196 unreducedVarID_.resize(vars.size());
00197
00198
00199
00200
00201 for (int b=0; b<vars.size(); b++)
00202 {
00203 Tabs tab2;
00204 unreducedVarID_[b].resize(vars[b].size());
00205 int k=0;
00206 for (int i=0; i<vars[b].size(); i++)
00207 {
00208 const FuncElementBase* t
00209 = dynamic_cast<const FuncElementBase*>(vars[b][i].ptr().get());
00210 int fid = t->fid().dofID();
00211 if (varFuncSet_.contains(fid)) continue;
00212 varFuncData_[b].append(getSharedFunctionData(t));
00213 varFuncSet_.put(fid);
00214 varIDToBlockMap_.put(fid, b);
00215 varIDToReducedIDMap_[b].put(fid, k);
00216 unreducedVarID_[b][k] = fid;
00217 k++;
00218 }
00219 SUNDANCE_MSG2(verb, tab2 << "block=" << b << " var functions are "
00220 << unreducedVarID_[b]);
00221 }
00222
00223
00224 unkIDToReducedIDMap_.resize(unks.size());
00225 unreducedUnkID_.resize(unks.size());
00226 for (int b=0; b<unks.size(); b++)
00227 {
00228 Tabs tab2;
00229 unreducedUnkID_[b].resize(unks[b].size());
00230 int k=0;
00231 for (int i=0; i<unks[b].size(); i++)
00232 {
00233 const UnknownFuncElement* u
00234 = dynamic_cast<const UnknownFuncElement*>(unks[b][i].ptr().get());
00235 TEUCHOS_TEST_FOR_EXCEPTION(u==0, std::runtime_error,
00236 "EquationSet ctor input unk function "
00237 << unks[b][i]
00238 << " does not appear to be a unk function");
00239 int fid = u->fid().dofID();
00240 if (unkFuncSet_.contains(fid)) continue;
00241 unkFuncData_[b].append(getSharedFunctionData(u));
00242 unkFuncSet_.put(fid);
00243 unkIDToBlockMap_.put(fid, b);
00244 unkIDToReducedIDMap_[b].put(fid, k);
00245 unreducedUnkID_[b][k] = fid;
00246 k++;
00247 }
00248 SUNDANCE_MSG2(verb, tab2 << "block=" << b << " unk functions are "
00249 << unreducedUnkID_[b]);
00250 }
00251
00252
00253
00254
00255
00256 unreducedUnkParamID_.resize(unkParams.size());
00257 for (int i=0; i<unkParams.size(); i++)
00258 {
00259 const UnknownParameterElement* u
00260 = dynamic_cast<const UnknownParameterElement*>(unkParams[i].ptr().get());
00261 TEUCHOS_TEST_FOR_EXCEPTION(u==0, std::runtime_error,
00262 "EquationSet ctor input unk parameter "
00263 << unkParams[i]
00264 << " does not appear to be a unk parameter");
00265 int fid = u->fid().dofID();
00266 unkParamSet_.put(fid);
00267 unkParamIDToReducedUnkParamIDMap_.put(fid, i);
00268 unreducedUnkParamID_[i] = fid;
00269 }
00270 SUNDANCE_MSG2(verb, tab1 << "unk parameters are "
00271 << unreducedUnkParamID_);
00272
00273
00274
00275
00276 unreducedFixedParamID_.resize(fixedParams.size());
00277 for (int i=0; i<fixedParams.size(); i++)
00278 {
00279 const UnknownParameterElement* u
00280 = dynamic_cast<const UnknownParameterElement*>(fixedParams[i].ptr().get());
00281 TEUCHOS_TEST_FOR_EXCEPTION(u==0, std::runtime_error,
00282 "EquationSet ctor input fixed parameter "
00283 << fixedParams[i]
00284 << " does not appear to be a fixed parameter");
00285 int fid = u->fid().dofID();
00286 fixedParamSet_.put(fid);
00287 fixedParamIDToReducedFixedParamIDMap_.put(fid, i);
00288 unreducedFixedParamID_[i] = fid;
00289 }
00290 SUNDANCE_MSG2(verb, tab1 << "fixed parameters are "
00291 << unreducedFixedParamID_);
00292
00293 Set<OrderedHandle<CellFilterStub> > regionSet;
00294
00295
00296 SUNDANCE_MSG1(verb, tab1 << "processing integral terms");
00297 {
00298 Tabs tab2;
00299 SUNDANCE_MSG3(verb, tab2 << integralSum_->rqcToExprMap());
00300 }
00301 for (Sundance::Map<RegionQuadCombo, Expr>::const_iterator
00302 r=integralSum_->rqcToExprMap().begin();
00303 r!=integralSum_->rqcToExprMap().end(); r++)
00304 {
00305 Tabs tab15;
00306 Tabs tab2;
00307 RegionQuadCombo rqc = r->first;
00308 int rqcVerb = verb;
00309 if (rqc.watch().isActive())
00310 {
00311 rqcVerb=rqc.watch().param("equation set setup");
00312 SUNDANCE_MSG1(max(verb, rqcVerb), tab15 << "processing RQC = " << rqc);
00313 }
00314
00315 Expr term = r->second;
00316
00317 SUNDANCE_MSG2(rqcVerb, tab2 << "expr = " << term);
00318 OrderedHandle<CellFilterStub> reg = rqc.domain();
00319
00320 if (!regionSet.contains(reg))
00321 {
00322 regionSet.put(reg);
00323 Set<int> vf = integralSum_->funcsOnRegion(reg, varFuncSet_);
00324 Set<int> uf = integralSum_->funcsOnRegion(reg, unkFuncSet_);
00325 SUNDANCE_MSG2(rqcVerb, tab2 << "vf = " << vf);
00326 SUNDANCE_MSG2(rqcVerb, tab2 << "uf = " << uf);
00327 varsOnRegions_.put(reg, vf);
00328 unksOnRegions_.put(reg, uf);
00329 }
00330 else
00331 {
00332 Set<int>& t = varsOnRegions_.get(reg);
00333 t.merge(integralSum_->funcsOnRegion(reg, varFuncSet_));
00334 Set<int>& u = unksOnRegions_.get(reg);
00335 u.merge(integralSum_->funcsOnRegion(reg, unkFuncSet_));
00336 }
00337 SUNDANCE_MSG2(rqcVerb, tab2 << "varsOnRegion = " << varsOnRegions_.get(reg));
00338 SUNDANCE_MSG2(rqcVerb, tab2 << "unksOnRegion = " << unksOnRegions_.get(reg));
00339 }
00340
00341
00342
00343 if (hasBCs)
00344 {
00345
00346
00347 for (Sundance::Map<RegionQuadCombo, Expr>::const_iterator
00348 r=bcSum_->rqcToExprMap().begin();
00349 r!=bcSum_->rqcToExprMap().end(); r++)
00350 {
00351 Tabs tab15;
00352 RegionQuadCombo rqc = r->first;
00353 int rqcVerb = verb;
00354 if (rqc.watch().isActive())
00355 {
00356 rqcVerb=rqc.watch().param("equation set setup");
00357 SUNDANCE_MSG1(max(verb, rqcVerb), tab15 << "processing RQC = " << rqc);
00358 }
00359
00360 Expr term = r->second;
00361 OrderedHandle<CellFilterStub> reg = rqc.domain();
00362
00363 if (!regionSet.contains(reg))
00364 {
00365 regionSet.put(reg);
00366 varsOnRegions_.put(reg, bcSum_->funcsOnRegion(reg, varFuncSet_));
00367 unksOnRegions_.put(reg, bcSum_->funcsOnRegion(reg, unkFuncSet_));
00368 bcVarsOnRegions_.put(reg, bcSum_->funcsOnRegion(reg, varFuncSet_));
00369 bcUnksOnRegions_.put(reg, bcSum_->funcsOnRegion(reg, unkFuncSet_));
00370 }
00371 else
00372 {
00373 if (!bcVarsOnRegions_.containsKey(reg))
00374 {
00375 bcVarsOnRegions_.put(reg, bcSum_->funcsOnRegion(reg, varFuncSet_));
00376 }
00377 if (!bcUnksOnRegions_.containsKey(reg))
00378 {
00379 bcUnksOnRegions_.put(reg, bcSum_->funcsOnRegion(reg, unkFuncSet_));
00380 }
00381 Set<int>& t = varsOnRegions_.get(reg);
00382 t.merge(bcSum_->funcsOnRegion(reg, varFuncSet_));
00383 Set<int>& u = unksOnRegions_.get(reg);
00384 u.merge(bcSum_->funcsOnRegion(reg, unkFuncSet_));
00385 }
00386 }
00387 }
00388
00389
00390 regions_ = regionSet.elements();
00391
00392 reducedVarsOnRegions_.resize(regions_.size());
00393 reducedUnksOnRegions_.resize(regions_.size());
00394 for (int r=0; r<regions_.size(); r++)
00395 {
00396 Tabs tab1;
00397 regionToIndexMap_.put(regions_[r], r);
00398 reducedVarsOnRegions_[r].resize(numVarBlocks());
00399 reducedUnksOnRegions_[r].resize(numUnkBlocks());
00400 OrderedHandle<CellFilterStub> cf = regions_[r];
00401 const Set<int>& v = this->varsOnRegion(r);
00402 const Set<int>& u = this->unksOnRegion(r);
00403 const Set<int>& bv = this->varsOnRegion(r);
00404 const Set<int>& bu = this->unksOnRegion(r);
00405 Set<int> vf = v;
00406 Set<int> uf = u;
00407 vf.merge(bv);
00408 uf.merge(bu);
00409 for (Set<int>::const_iterator i=vf.begin(); i!=vf.end(); i++)
00410 {
00411 int fid = *i;
00412 SUNDANCE_MSG2(verb, tab1 << "adding VF=" << fid << " to testToRegionMap");
00413 if (testToRegionsMap_.containsKey(fid))
00414 {
00415 testToRegionsMap_[fid].put(cf);
00416 }
00417 else
00418 {
00419 Set<OrderedHandle<CellFilterStub> > s;
00420 s.put(cf);
00421 testToRegionsMap_.put(fid, s);
00422 }
00423 int rv = reducedVarID(fid);
00424 int br = blockForVarID(fid);
00425 reducedVarsOnRegions_[r][br].put(rv);
00426 }
00427 for (Set<int>::const_iterator i=uf.begin(); i!=uf.end(); i++)
00428 {
00429 int fid = *i;
00430 if (unkToRegionsMap_.containsKey(fid))
00431 {
00432 unkToRegionsMap_[fid].put(cf);
00433 }
00434 else
00435 {
00436 Set<OrderedHandle<CellFilterStub> > s;
00437 s.put(cf);
00438 unkToRegionsMap_.put(fid, s);
00439 }
00440 int ru = reducedUnkID(fid);
00441 int bc = blockForUnkID(fid);
00442 reducedUnksOnRegions_[r][bc].put(ru);
00443 }
00444 }
00445 }
00446
00447
00448
00449
00450 Array<Expr> FunctionSupportResolver::flattenSpectral(const Array<Expr>& expr) const
00451 {
00452 Array<Expr> rtn(expr.size());
00453 for (int i=0; i<expr.size(); i++)
00454 {
00455 const Expr& e = expr[i];
00456 rtn[i] = flattenSpectral(e);
00457 }
00458 return rtn;
00459 }
00460
00461 Expr FunctionSupportResolver::flattenSpectral(const Expr& expr) const
00462 {
00463 Array<Expr> rtn(expr.size());
00464 for (int i=0; i<expr.size(); i++)
00465 {
00466 if (expr[i].size() == 1)
00467 {
00468 const SpectralExpr* se
00469 = dynamic_cast<const SpectralExpr*>(expr[i][0].ptr().get());
00470 if (se != 0)
00471 {
00472 int nt = se->getSpectralBasis().nterms();
00473 Array<Expr> e(nt);
00474 for (int j=0; j<nt; j++)
00475 {
00476 e[j] = se->getCoeff(j);
00477 }
00478 rtn[i] = new ListExpr(e);
00479 }
00480 else
00481 {
00482 rtn[i] = expr[i];
00483 }
00484 }
00485 else
00486 {
00487 rtn[i] = flattenSpectral(expr[i]);
00488 }
00489 }
00490 Expr r = new ListExpr(rtn);
00491 return r.flatten();
00492
00493 }
00494
00495 int FunctionSupportResolver::reducedVarID(int varID) const
00496 {
00497 TEUCHOS_TEST_FOR_EXCEPTION(!hasVarID(varID), std::runtime_error,
00498 "varID " << varID << " not found in equation set");
00499
00500 int b = blockForVarID(varID);
00501 return varIDToReducedIDMap_[b].get(varID);
00502 }
00503
00504 int FunctionSupportResolver::reducedUnkID(int unkID) const
00505 {
00506 TEUCHOS_TEST_FOR_EXCEPTION(!hasUnkID(unkID), std::runtime_error,
00507 "unkID " << unkID << " not found in equation set");
00508
00509 int b = blockForUnkID(unkID);
00510 return unkIDToReducedIDMap_[b].get(unkID);
00511 }
00512
00513
00514 int FunctionSupportResolver::reducedUnkParamID(int unkID) const
00515 {
00516 TEUCHOS_TEST_FOR_EXCEPTION(!hasUnkParamID(unkID), std::runtime_error,
00517 "unkParamID " << unkID << " not found in equation set");
00518
00519 return unkParamIDToReducedUnkParamIDMap_.get(unkID);
00520 }
00521
00522 int FunctionSupportResolver::reducedFixedParamID(int parID) const
00523 {
00524 TEUCHOS_TEST_FOR_EXCEPTION(!hasFixedParamID(parID), std::runtime_error,
00525 "fixedParamID " << parID << " not found in equation set");
00526
00527 return fixedParamIDToReducedFixedParamIDMap_.get(parID);
00528 }
00529
00530 int FunctionSupportResolver::blockForVarID(int varID) const
00531 {
00532 TEUCHOS_TEST_FOR_EXCEPTION(!varIDToBlockMap_.containsKey(varID), std::runtime_error,
00533 "key " << varID << " not found in map "
00534 << varIDToBlockMap_);
00535 return varIDToBlockMap_.get(varID);
00536 }
00537
00538 int FunctionSupportResolver::blockForUnkID(int unkID) const
00539 {
00540 TEUCHOS_TEST_FOR_EXCEPTION(!unkIDToBlockMap_.containsKey(unkID), std::runtime_error,
00541 "key " << unkID << " not found in map "
00542 << unkIDToBlockMap_);
00543 return unkIDToBlockMap_.get(unkID);
00544 }
00545
00546 const Set<OrderedHandle<CellFilterStub> >& FunctionSupportResolver::regionsForTestFunc(int testID) const
00547 {
00548 TEUCHOS_TEST_FOR_EXCEPTION(!testToRegionsMap_.containsKey(testID), std::runtime_error,
00549 "key " << testID << " not found in map "
00550 << testToRegionsMap_);
00551 return testToRegionsMap_.get(testID);
00552 }
00553
00554 const Set<OrderedHandle<CellFilterStub> >& FunctionSupportResolver::regionsForUnkFunc(int unkID) const
00555 {
00556 TEUCHOS_TEST_FOR_EXCEPTION(!unkToRegionsMap_.containsKey(unkID), std::runtime_error,
00557 "key " << unkID << " not found in map "
00558 << testToRegionsMap_);
00559 return unkToRegionsMap_.get(unkID);
00560 }
00561
00562 int FunctionSupportResolver::indexForRegion(const OrderedHandle<CellFilterStub>& region) const
00563 {
00564 TEUCHOS_TEST_FOR_EXCEPTION(!regionToIndexMap_.containsKey(region), std::runtime_error,
00565 "key " << region << " not found in map "
00566 << regionToIndexMap_);
00567 return regionToIndexMap_.get(region);
00568 }
00569
00570 bool FunctionSupportResolver::hasBCs() const
00571 {
00572 return bcSum_ != 0;
00573 }
00574
00575
00576 namespace Sundance
00577 {
00578
00579 RCP<const CommonFuncDataStub> getSharedFunctionData(const FuncElementBase* f)
00580 {
00581 const UnknownFuncElement* u = dynamic_cast<const UnknownFuncElement*>(f);
00582 const TestFuncElement* t = dynamic_cast<const TestFuncElement*>(f);
00583
00584 if (u != 0)
00585 {
00586 return rcp_dynamic_cast<const CommonFuncDataStub>(u->commonData());
00587 }
00588 if (t != 0)
00589 {
00590 return rcp_dynamic_cast<const CommonFuncDataStub>(t->commonData());
00591 }
00592
00593 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error,
00594 "unrecognized function type: " << typeid(*f).name());
00595 return u->commonData();
00596 }
00597 }
00598