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 "SundanceDOFMapBuilder.hpp"
00032 #include "SundanceOut.hpp"
00033 #include "PlayaTabs.hpp"
00034 #include "SundanceBasisFamily.hpp"
00035 #include "SundanceLagrange.hpp"
00036 #include "SundanceEdgeLocalizedBasis.hpp"
00037 #include "SundanceMixedDOFMap.hpp"
00038 #include "SundanceMixedDOFMapHN.hpp"
00039 #include "SundanceNodalDOFMap.hpp"
00040 #include "SundanceSubmaximalNodalDOFMap.hpp"
00041 #include "SundanceNodalDOFMapHN.hpp"
00042 #include "SundancePartialElementDOFMap.hpp"
00043 #include "SundanceMaximalCellFilter.hpp"
00044 #include "SundanceInhomogeneousNodalDOFMap.hpp"
00045 #include "SundanceInhomogeneousEdgeLocalizedDOFMap.hpp"
00046 #include "SundanceInhomogeneousDOFMapHN.hpp"
00047 #include "SundanceCellFilter.hpp"
00048 #include "SundanceDimensionalCellFilter.hpp"
00049 #include "SundanceCellSet.hpp"
00050 #include "SundanceCFMeshPair.hpp"
00051 #include "Teuchos_Time.hpp"
00052 #include "Teuchos_TimeMonitor.hpp"
00053
00054
00055 using namespace Sundance;
00056 using namespace Teuchos;
00057 using Playa::MPIComm;
00058 using Playa::MPIDataType;
00059 using Playa::MPIOp;
00060
00061 static Time& DOFBuilderCtorTimer()
00062 {
00063 static RCP<Time> rtn
00064 = TimeMonitor::getNewTimer("DOF map building");
00065 return *rtn;
00066 }
00067
00068 static Time& cellFilterReductionTimer()
00069 {
00070 static RCP<Time> rtn
00071 = TimeMonitor::getNewTimer("cell filter reduction");
00072 return *rtn;
00073 }
00074
00075 static Time& findFuncDomainTimer()
00076 {
00077 static RCP<Time> rtn
00078 = TimeMonitor::getNewTimer("finding func domains");
00079 return *rtn;
00080 }
00081
00082 DOFMapBuilder::DOFMapBuilder(const Mesh& mesh,
00083 const RCP<FunctionSupportResolver>& fsr, bool findBCCols,
00084 int setupVerb)
00085 : verb_(setupVerb),
00086 mesh_(mesh),
00087 fsr_(fsr),
00088 rowMap_(),
00089 colMap_(),
00090 isBCRow_(),
00091 isBCCol_(),
00092 remoteBCCols_()
00093 {
00094 init(findBCCols);
00095 }
00096
00097 DOFMapBuilder::DOFMapBuilder(int setupVerb)
00098 : verb_(setupVerb),
00099 mesh_(),
00100 fsr_(),
00101 rowMap_(),
00102 colMap_(),
00103 isBCRow_(),
00104 isBCCol_(),
00105 remoteBCCols_()
00106 {}
00107
00108 RCP<DOFMapBase> DOFMapBuilder::makeMap(const Mesh& mesh,
00109 const Array<RCP<BasisDOFTopologyBase> >& basis,
00110 const Array<Set<CellFilter> >& filters)
00111 {
00112 verb_=0;
00113 TimeMonitor timer(DOFBuilderCtorTimer());
00114 SUNDANCE_MSG1(verb_, "in DOFMapBuilder::makeMap()");
00115 for (int i=0; i<basis.size(); i++)
00116 {
00117 SUNDANCE_MSG2(verb_, "i=" << i
00118 << " basis=" << basis[i]->description()
00119 << " filters=" << filters[i]);
00120 }
00121
00122 RCP<DOFMapBase> rtn;
00123
00124 if (allowNodalMap() && hasOmnipresentNodalMap(basis, mesh, filters))
00125 {
00126 SUNDANCE_MSG2(verb_, "creating omnipresent nodal map");
00127 CellFilter maxCells = getMaxCellFilter(filters);
00128
00129 if (mesh.allowsHangingHodes()){
00130 rtn = rcp(new NodalDOFMapHN(mesh, basis.size(), maxCells, verb_));
00131 }
00132 else {
00133 rtn = rcp(new NodalDOFMap(mesh, basis.size(), maxCells, verb_));
00134 }
00135 }
00136 else if (hasNodalBasis(basis) && filtersAreZeroDimensional(mesh, filters))
00137 {
00138 SUNDANCE_MSG2(verb_, "creating submaximal nodal map");
00139 TEUCHOS_TEST_FOR_EXCEPT(filters.size() != 1);
00140 TEUCHOS_TEST_FOR_EXCEPT(filters[0].size() != 1);
00141 rtn = rcp(new SubmaximalNodalDOFMap(mesh, *filters[0].begin(), basis.size(), verb_));
00142 }
00143 else if (hasCellBasis(basis) && hasCommonDomain(filters))
00144 {
00145 TEUCHOS_TEST_FOR_EXCEPTION(filters[0].size() != 1, std::runtime_error,
00146 "only a single domain expected in construction of an element "
00147 "DOF map");
00148 rtn = rcp(new PartialElementDOFMap(mesh, *filters[0].begin(), basis.size(), verb_));
00149 }
00150 else if (allFuncsAreOmnipresent(mesh, filters))
00151 {
00152 SUNDANCE_MSG2(verb_, "creating omnipresent mixed map");
00153 CellFilter maxCells = getMaxCellFilter(filters);
00154 if (mesh.allowsHangingHodes()){
00155 rtn = rcp(new MixedDOFMapHN(mesh, basis, maxCells, verb_));
00156 }else
00157 {
00158 rtn = rcp(new MixedDOFMap(mesh, basis, maxCells, verb_));
00159 }
00160 }
00161 else if ( hasNodalBasis(basis) && (!mesh.allowsHangingHodes()) )
00162 {
00163 SUNDANCE_MSG2(verb_, "creating inhomogeneous nodal map");
00164 Sundance::Map<CellFilter, Set<int> > fmap = domainToFuncSetMap(filters);
00165 Sundance::Map<CellFilter, Sundance::Map<Set<int>, CellSet> > inputChildren;
00166
00167 Array<Sundance::Map<Set<int>, CellFilter> > disjoint
00168 = DOFMapBuilder::funcDomains(mesh, fmap, inputChildren);
00169
00170 rtn = rcp(new InhomogeneousNodalDOFMap(mesh, disjoint, verb_));
00171 }
00172 else if ( hasEdgeLocalizedBasis(basis) && (!mesh.allowsHangingHodes()) )
00173 {
00174 SUNDANCE_MSG2(verb_, "creating inhomogeneous edge-localized map");
00175
00176 Sundance::Map<CellFilter, Set<int> > fmap = domainToFuncSetMap(filters);
00177
00178 Sundance::Map<CellFilter, Sundance::Map<Set<int>, CellSet> > inputChildren;
00179 Array<Sundance::Map<Set<int>, CellFilter> > disjoint
00180 = DOFMapBuilder::funcDomains(mesh, fmap, inputChildren);
00181
00182 rtn = rcp(new InhomogeneousEdgeLocalizedDOFMap(mesh, disjoint, verb_));
00183 }
00184 else
00185 {
00186 SUNDANCE_MSG2(verb_, "creating inhomogeneous HN map (the last possibility)");
00187 Sundance::Map<CellFilter, Set<int> > fmap = domainToFuncSetMap(filters);
00188 Sundance::Map<CellFilter, Sundance::Map<Set<int>, CellSet> > inputChildren;
00189
00190 Array<Sundance::Map<Set<int>, CellFilter> > disjoint
00191 = DOFMapBuilder::funcDomains(mesh, fmap, inputChildren);
00192
00193 rtn = rcp(new InhomogeneousDOFMapHN(mesh, basis , disjoint, verb_));
00194 }
00195
00196 if (verb_ > 0)
00197 {
00198 Out::os() << "done building DOF map" << std::endl;
00199 if (verb_ > 1) Out::os() << "num DOFs" << rtn->numDOFs() << std::endl;
00200 if (verb_ > 1) Out::os() << "num local DOFs"
00201 << rtn->numLocalDOFs() << std::endl;
00202 if (verb_ > 4) rtn->print(Out::os());
00203 }
00204 return rtn;
00205 }
00206
00207
00208 Sundance::Map<CellFilter, Set<int> > DOFMapBuilder::domainToFuncSetMap(const Array<Set<CellFilter> >& filters) const
00209 {
00210 SUNDANCE_MSG2(verb_, "in DOFMapBuilder::domainToFuncSetMap()");
00211 Map<CellFilter, Set<int> > rtn;
00212 for (int i=0; i<filters.size(); i++)
00213 {
00214 const Set<CellFilter>& s = filters[i];
00215 for (Set<CellFilter>::const_iterator j=s.begin(); j!=s.end(); j++)
00216 {
00217 const CellFilter& cf = *j;
00218 if (rtn.containsKey(cf))
00219 {
00220 rtn[cf].put(i);
00221 }
00222 else
00223 {
00224
00225 rtn.put(cf, makeSet((int) i));
00226 }
00227 }
00228 }
00229 for (Map<CellFilter, Set<int> >::const_iterator
00230 i=rtn.begin(); i!=rtn.end(); i++)
00231 {
00232 SUNDANCE_MSG2(verb_, "subdomain=" << i->first << ", functions="
00233 << i->second);
00234 }
00235 return rtn;
00236 }
00237
00238
00239 void DOFMapBuilder
00240 ::getSubdomainUnkFuncMatches(const FunctionSupportResolver& fsr,
00241 Array<Sundance::Map<CellFilter, Set<int> > >& fmap) const
00242 {
00243 fmap.resize(fsr.numUnkBlocks());
00244
00245 for (int r=0; r<fsr.numRegions(); r++)
00246 {
00247 CellFilter subreg = fsr.region(r);
00248 Set<int> funcs = fsr.unksOnRegion(r).setUnion(fsr.bcUnksOnRegion(r));
00249 for (Set<int>::const_iterator i=funcs.begin(); i!=funcs.end(); i++)
00250 {
00251 int block = fsr.blockForUnkID(*i);
00252 if (fmap[block].containsKey(subreg))
00253 {
00254 fmap[block][subreg].put(*i);
00255 }
00256 else
00257 {
00258 fmap[block].put(subreg, makeSet(*i));
00259 }
00260 }
00261 }
00262 }
00263
00264 void DOFMapBuilder
00265 ::getSubdomainVarFuncMatches(const FunctionSupportResolver& fsr,
00266 Array<Sundance::Map<CellFilter, Set<int> > >& fmap) const
00267 {
00268 fmap.resize(fsr.numVarBlocks());
00269
00270 for (int r=0; r<fsr.numRegions(); r++)
00271 {
00272 CellFilter subreg = fsr.region(r);
00273 Set<int> funcs = fsr.varsOnRegion(r).setUnion(fsr.bcVarsOnRegion(r));
00274 for (Set<int>::const_iterator i=funcs.begin(); i!=funcs.end(); i++)
00275 {
00276 int block = fsr.blockForVarID(*i);
00277 if (fmap[block].containsKey(subreg))
00278 {
00279 fmap[block][subreg].put(*i);
00280 }
00281 else
00282 {
00283 fmap[block].put(subreg, makeSet(*i));
00284 }
00285 }
00286 }
00287 }
00288
00289 Array<Sundance::Map<Set<int>, CellFilter> > DOFMapBuilder
00290 ::funcDomains(const Mesh& mesh,
00291 const Sundance::Map<CellFilter, Set<int> >& fmap,
00292 Sundance::Map<CellFilter, Sundance::Map<Set<int>, CellSet> >& inputToChildrenMap) const
00293 {
00294 TimeMonitor timer(findFuncDomainTimer());
00295 Array<Array<CellFilter> > filters(mesh.spatialDim()+1);
00296 Array<Array<Set<int> > > funcs(mesh.spatialDim()+1);
00297
00298 for (Sundance::Map<CellFilter, Set<int> >::const_iterator
00299 i=fmap.begin(); i!=fmap.end(); i++)
00300 {
00301 int d = i->first.dimension(mesh);
00302 filters[d].append(i->first);
00303 funcs[d].append(i->second);
00304 }
00305 Array<Array<CFMeshPair> > tmp(mesh.spatialDim()+1);
00306 for (int d=0; d<tmp.size(); d++)
00307 {
00308 if (filters[d].size() != 0)
00309 tmp[d] = findDisjointFilters(filters[d], funcs[d], mesh);
00310 }
00311
00312 for (int d=0; d<tmp.size(); d++)
00313 {
00314 for (int r=0; r<tmp[d].size(); r++)
00315 {
00316 for (int p=0; p<filters[d].size(); p++)
00317 {
00318 if (tmp[d][r].filter().isSubsetOf(filters[d][p], mesh))
00319 {
00320 if (inputToChildrenMap.containsKey(filters[d][p]))
00321 {
00322 Sundance::Map<Set<int>, CellSet>& m
00323 = inputToChildrenMap[filters[d][p]];
00324 if (m.containsKey(tmp[d][r].funcs()))
00325 {
00326 m.put(tmp[d][r].funcs(), m[tmp[d][r].funcs()].setUnion(tmp[d][r].cellSet()));
00327 }
00328 else
00329 {
00330 m.put(tmp[d][r].funcs(), tmp[d][r].cellSet());
00331 }
00332 }
00333 else
00334 {
00335 Sundance::Map<Set<int>, CellSet> m;
00336 m.put(tmp[d][r].funcs(), tmp[d][r].cellSet());
00337 inputToChildrenMap.put(filters[d][p], m);
00338 }
00339 }
00340 }
00341 }
00342 }
00343
00344 Array<Sundance::Map<Set<int>, CellFilter> > rtn(mesh.spatialDim()+1);
00345 for (int d=0; d<tmp.size(); d++)
00346 {
00347 if (tmp[d].size() == 0) continue;
00348 for (int i=0; i<tmp[d].size(); i++)
00349 {
00350 const Set<int>& f = tmp[d][i].funcs();
00351 const CellFilter& cf = tmp[d][i].filter();
00352 if (rtn[d].containsKey(f))
00353 {
00354 rtn[d].put(f, rtn[d][f] + cf);
00355 }
00356 else
00357 {
00358 rtn[d].put(f, cf);
00359 }
00360 }
00361 }
00362
00363 return rtn;
00364 }
00365
00366
00367 void DOFMapBuilder::init(bool findBCCols)
00368 {
00369 SUNDANCE_MSG1(verb_, "in DOFMapBuilder::init()");
00370 SUNDANCE_MSG2(verb_, "num var blocks=" << fsr_->numVarBlocks());
00371 SUNDANCE_MSG2(verb_, "num unk blocks=" << fsr_->numUnkBlocks());
00372
00373 rowMap_.resize(fsr_->numVarBlocks());
00374 colMap_.resize(fsr_->numUnkBlocks());
00375 isBCRow_.resize(fsr_->numVarBlocks());
00376 isBCCol_.resize(fsr_->numUnkBlocks());
00377
00378 Array<Array<RCP<BasisDOFTopologyBase> > > testBasis = testBasisTopologyArray();
00379 Array<Array<Set<CellFilter> > > testRegions = testCellFilters();
00380
00381 for (int br=0; br<fsr_->numVarBlocks(); br++)
00382 {
00383 SUNDANCE_MSG2(verb_, "making map for block row=" << br);
00384 rowMap_[br] = makeMap(mesh_, testBasis[br], testRegions[br]);
00385 SUNDANCE_MSG2(verb_, "marking BC rows for block row=" << br);
00386 markBCRows(br);
00387 }
00388
00389
00390 Array<Array<RCP<BasisDOFTopologyBase> > > unkBasis = unkBasisTopologyArray();
00391 Array<Array<Set<CellFilter> > > unkRegions = unkCellFilters();
00392
00393 for (int bc=0; bc<fsr_->numUnkBlocks(); bc++)
00394 {
00395 if (isSymmetric(bc))
00396 {
00397 colMap_[bc] = rowMap_[bc];
00398 }
00399 else
00400 {
00401 SUNDANCE_MSG2(verb_, "making map for block col=" << bc);
00402 colMap_[bc] = makeMap(mesh_, unkBasis[bc], unkRegions[bc]);
00403 }
00404 SUNDANCE_MSG2(verb_, "marking BC cols for block col=" << bc);
00405 if (findBCCols) markBCCols(bc);
00406 }
00407 }
00408
00409 void DOFMapBuilder::extractUnkSetsFromFSR(const FunctionSupportResolver& fsr,
00410 Array<Set<int> >& funcSets,
00411 Array<CellFilter>& regions) const
00412 {
00413 funcSets.resize(fsr.numRegions());
00414 regions.resize(fsr.numRegions());
00415 for (int r=0; r<fsr.numRegions(); r++)
00416 {
00417 regions[r] = fsr.region(r);
00418 funcSets[r] = fsr.unksOnRegion(r).setUnion(fsr.bcUnksOnRegion(r));
00419 }
00420 }
00421
00422 void DOFMapBuilder::extractVarSetsFromFSR(const FunctionSupportResolver& fsr,
00423 Array<Set<int> >& funcSets,
00424 Array<CellFilter>& regions) const
00425 {
00426 funcSets.resize(fsr.numRegions());
00427 regions.resize(fsr.numRegions());
00428 for (int r=0; r<fsr.numRegions(); r++)
00429 {
00430 regions[r] = fsr.region(r);
00431 funcSets[r] = fsr.varsOnRegion(r).setUnion(fsr.bcVarsOnRegion(r));
00432 }
00433 }
00434
00435 Sundance::Map<Set<int>, Set<CellFilter> >
00436 DOFMapBuilder::buildFuncSetToCFSetMap(const Array<Set<int> >& funcSets,
00437 const Array<CellFilter>& regions,
00438 const Mesh& mesh) const
00439 {
00440 Sundance::Map<Set<int>, Set<CellFilter> > tmp;
00441
00442 for (int r=0; r<regions.size(); r++)
00443 {
00444 const CellFilter& reg = regions[r];
00445 if (!tmp.containsKey(funcSets[r]))
00446 {
00447 tmp.put(funcSets[r], Set<CellFilter>());
00448 }
00449 tmp[funcSets[r]].put(reg);
00450 }
00451
00452
00453
00454 Sundance::Map<Set<int>, Set<CellFilter> > rtn=tmp;
00455
00456
00457
00458
00459
00460
00461
00462 return rtn;
00463 }
00464
00465 bool DOFMapBuilder::hasOmnipresentNodalMap(const Array<RCP<BasisDOFTopologyBase> >& basis,
00466 const Mesh& mesh,
00467 const Array<Set<CellFilter> >& filters) const
00468 {
00469 return hasNodalBasis(basis) && allFuncsAreOmnipresent(mesh, filters);
00470 }
00471
00472
00473
00474 bool DOFMapBuilder::hasCommonDomain(const Array<Set<CellFilter> >& filters) const
00475 {
00476 Set<CellFilter> first = filters[0];
00477 for (int i=1; i<filters.size(); i++)
00478 {
00479 if (! (filters[i] == first) ) return false;
00480 }
00481 return true;
00482 }
00483
00484 bool DOFMapBuilder::hasNodalBasis(const Array<RCP<BasisDOFTopologyBase> >& basis) const
00485 {
00486 for (int i=0; i<basis.size(); i++)
00487 {
00488 const Lagrange* lagr
00489 = dynamic_cast<const Lagrange*>(basis[i].get());
00490 if (lagr==0 || lagr->order()!=1) return false;
00491 }
00492 return true;
00493 }
00494
00495 bool DOFMapBuilder::hasEdgeLocalizedBasis(const Array<RCP<BasisDOFTopologyBase> >& basis) const
00496 {
00497 for (int i=0; i<basis.size(); i++)
00498 {
00499 const RCP<const EdgeLocalizedBasis> downcasted = rcp_dynamic_cast<EdgeLocalizedBasis>(basis[i]);
00500 if (is_null(downcasted)) return false;
00501 }
00502 return true;
00503 }
00504
00505 bool DOFMapBuilder::hasCellBasis(const Array<RCP<BasisDOFTopologyBase> >& basis) const
00506 {
00507 for (int i=0; i<basis.size(); i++)
00508 {
00509 const Lagrange* lagr
00510 = dynamic_cast<const Lagrange*>(basis[i].get());
00511 if (lagr==0 || lagr->order()!=0) return false;
00512 }
00513 return true;
00514 }
00515
00516 bool DOFMapBuilder::filtersAreZeroDimensional(const Mesh& mesh,
00517 const Array<Set<CellFilter> >& filters) const
00518 {
00519 for (int i=0; i<filters.size(); i++)
00520 {
00521 for (Set<CellFilter>::const_iterator
00522 j=filters[i].begin(); j!=filters[i].end(); j++)
00523 {
00524 if (j->dimension(mesh) != 0) return false;
00525 }
00526 }
00527 return true;
00528 }
00529
00530
00531 bool DOFMapBuilder::allFuncsAreOmnipresent(const Mesh& mesh,
00532 const Array<Set<CellFilter> >& filters) const
00533 {
00534 int rtn = 0;
00535
00536 int maxFilterDim = 0;
00537 Set<Set<CellFilter> > distinctSets;
00538 for (int i=0; i<filters.size(); i++)
00539 {
00540 for (Set<CellFilter>::const_iterator iter=filters[i].begin();
00541 iter != filters[i].end(); iter++)
00542 {
00543 int dim = iter->dimension(mesh);
00544 if (dim > maxFilterDim) maxFilterDim = dim;
00545 }
00546 distinctSets.put(filters[i]);
00547 }
00548
00549 for (Set<Set<CellFilter> >::const_iterator
00550 iter=distinctSets.begin(); iter != distinctSets.end(); iter++)
00551 {
00552 if (!isWholeDomain(mesh, maxFilterDim, *iter)) rtn = 1;
00553 }
00554
00555
00556 int omniPresent = rtn;
00557 mesh.comm().allReduce((void*) &omniPresent, (void*) &rtn, 1,
00558 MPIDataType::intType(), MPIOp::sumOp());
00559
00560
00561 return (rtn < 1);
00562 }
00563
00564 bool DOFMapBuilder::isWholeDomain(const Mesh& mesh,
00565 int maxFilterDim,
00566 const Set<CellFilter>& filters) const
00567 {
00568 CellFilter allMax;
00569 if (maxFilterDim==mesh.spatialDim()) allMax = new MaximalCellFilter();
00570 else allMax = new DimensionalCellFilter(maxFilterDim);
00571
00572 CellSet remainder = allMax.getCells(mesh);
00573
00574 for (Set<CellFilter>::const_iterator
00575 i=filters.begin(); i!=filters.end(); i++)
00576 {
00577 const CellFilter& cf = *i;
00578 if (maxFilterDim==mesh.spatialDim())
00579 {
00580 if (0 != dynamic_cast<const MaximalCellFilter*>(cf.ptr().get()))
00581 {
00582 return true;
00583 }
00584 }
00585 else
00586 {
00587 const DimensionalCellFilter* dcf
00588 = dynamic_cast<const DimensionalCellFilter*>(cf.ptr().get());
00589 if (0 != dcf && dcf->dimension(mesh) == maxFilterDim)
00590 {
00591 return true;
00592 }
00593 }
00594 if (cf.dimension(mesh) != maxFilterDim) continue;
00595 CellSet cells = cf.getCells(mesh);
00596 SUNDANCE_MSG2(verb_, "found " << cells.numCells() << " in CF " << cf);
00597 remainder = remainder.setDifference(cells);
00598 if (remainder.begin() == remainder.end()) return true;
00599 }
00600
00601 SUNDANCE_MSG2(verb_, "num remaining cells: " << remainder.numCells());
00602
00603 return false;
00604 }
00605
00606
00607
00608 CellFilter DOFMapBuilder::getMaxCellFilter(const Array<Set<CellFilter> >& filters) const
00609 {
00610 for (int i=0; i<filters.size(); i++)
00611 {
00612 const Set<CellFilter>& cfs = filters[i];
00613 if (cfs.size() != 1) continue;
00614 const CellFilter& cf = *cfs.begin();
00615 if (0 != dynamic_cast<const MaximalCellFilter*>(cf.ptr().get()))
00616 return cf;
00617 }
00618
00619 return new MaximalCellFilter();
00620 }
00621
00622
00623
00624 Array<Array<Set<CellFilter> > > DOFMapBuilder::testCellFilters() const
00625 {
00626 Array<Array<Set<CellFilter> > > rtn(fsr_->numVarBlocks());
00627
00628 for (int b=0; b<fsr_->numVarBlocks(); b++)
00629 {
00630 for (int i=0; i<fsr_->numVarIDs(b); i++)
00631 {
00632 int testID = fsr_->unreducedVarID(b, i);
00633 Set<CellFilter> s;
00634 const Set<OrderedHandle<CellFilterStub> >& cfs
00635 = fsr_->regionsForTestFunc(testID);
00636 for (Set<OrderedHandle<CellFilterStub> >::const_iterator
00637 j=cfs.begin(); j!=cfs.end(); j++)
00638 {
00639 RCP<CellFilterBase> cfb
00640 = rcp_dynamic_cast<CellFilterBase>(j->ptr());
00641 TEUCHOS_TEST_FOR_EXCEPT(cfb.get()==0);
00642 CellFilter cf = j->ptr();
00643 s.put(cf);
00644 }
00645
00646
00647 rtn[b].append(s);
00648 }
00649 }
00650 return rtn;
00651 }
00652
00653 Array<Array<Set<CellFilter> > > DOFMapBuilder::unkCellFilters() const
00654 {
00655 Array<Array<Set<CellFilter> > > rtn(fsr_->numUnkBlocks());
00656
00657 for (int b=0; b<fsr_->numUnkBlocks(); b++)
00658 {
00659 for (int i=0; i<fsr_->numUnkIDs(b); i++)
00660 {
00661 int unkID = fsr_->unreducedUnkID(b, i);
00662 Set<CellFilter> s;
00663 const Set<OrderedHandle<CellFilterStub> >& cfs
00664 = fsr_->regionsForUnkFunc(unkID);
00665 for (Set<OrderedHandle<CellFilterStub> >::const_iterator
00666 j=cfs.begin(); j!=cfs.end(); j++)
00667 {
00668 RCP<CellFilterBase> cfb
00669 = rcp_dynamic_cast<CellFilterBase>(j->ptr());
00670 TEUCHOS_TEST_FOR_EXCEPT(cfb.get()==0);
00671 CellFilter cf = j->ptr();
00672 s.put(cf);
00673 }
00674
00675
00676 rtn[b].append(s);
00677 }
00678 }
00679 return rtn;
00680 }
00681
00682 Array<Array<RCP<BasisDOFTopologyBase> > > DOFMapBuilder::testBasisTopologyArray() const
00683 {
00684 Array<Array<RCP<BasisDOFTopologyBase> > > rtn(fsr_->numVarBlocks());
00685 for (int b=0; b<fsr_->numVarBlocks(); b++)
00686 {
00687 for (int i=0; i<fsr_->numVarIDs(b); i++)
00688 {
00689 rtn[b].append(BasisFamily::getBasisTopology(fsr_->varFuncData(b, i)));
00690 }
00691 }
00692 return rtn;
00693 }
00694
00695 Array<Array<RCP<BasisDOFTopologyBase> > > DOFMapBuilder::unkBasisTopologyArray() const
00696 {
00697 Array<Array<RCP<BasisDOFTopologyBase> > > rtn(fsr_->numUnkBlocks());
00698 for (int b=0; b<fsr_->numUnkBlocks(); b++)
00699 {
00700 for (int i=0; i<fsr_->numUnkIDs(b); i++)
00701 {
00702 rtn[b].append(BasisFamily::getBasisTopology(fsr_->unkFuncData(b, i)));
00703 }
00704 }
00705 return rtn;
00706 }
00707
00708
00709 Set<CellFilter> DOFMapBuilder
00710 ::reduceCellFilters(const Mesh& mesh,
00711 const Set<CellFilter>& inputSet) const
00712 {
00713 TimeMonitor timer(cellFilterReductionTimer());
00714 Set<CellFilter> rtn;
00715
00716 CellFilter m = new MaximalCellFilter();
00717 if (inputSet.contains(m))
00718 {
00719 rtn.put(m);
00720 return rtn;
00721 }
00722
00723
00724
00725 CellFilter myMaxFilters;
00726 for (Set<CellFilter>::const_iterator
00727 i=inputSet.begin(); i!=inputSet.end(); i++)
00728 {
00729 CellFilter f = *i;
00730 if (f.dimension(mesh) != mesh.spatialDim()) continue;
00731 myMaxFilters = myMaxFilters + f;
00732 }
00733
00734 CellSet myMax = myMaxFilters.getCells(mesh);
00735
00736
00737
00738
00739
00740
00741 CellSet allMax = m.getCells(mesh);
00742 CellSet diff = allMax.setDifference(myMax);
00743
00744
00745 if (diff.begin() == diff.end())
00746 {
00747 rtn.put(m);
00748 return rtn;
00749 }
00750
00751
00752
00753 if (myMax.begin() != myMax.end()) rtn.put(myMaxFilters);
00754
00755
00756
00757
00758
00759 for (Set<CellFilter>::const_iterator
00760 i=inputSet.begin(); i!=inputSet.end(); i++)
00761 {
00762 CellFilter f = *i;
00763 if (f.dimension(mesh) == mesh.spatialDim()) continue;
00764 CellSet s = f.getCells(mesh);
00765 if (s.areFacetsOf(myMax)) continue;
00766
00767
00768
00769
00770 rtn.put(f);
00771 }
00772 return rtn;
00773 }
00774
00775
00776 bool DOFMapBuilder::isSymmetric(int b) const
00777 {
00778 if ((int)fsr_->numVarBlocks() < b || (int)fsr_->numUnkBlocks() < b) return false;
00779
00780 if (fsr_->numVarIDs(b) != fsr_->numUnkIDs(b)) return false;
00781
00782 for (int i=0; i<fsr_->numVarIDs(b); i++)
00783 {
00784 BasisFamily basis1 = BasisFamily::getBasis(fsr_->varFuncData(b,i));
00785 BasisFamily basis2 = BasisFamily::getBasis(fsr_->unkFuncData(b,i));
00786 if (!(basis1 == basis2)) return false;
00787 }
00788 return true;
00789 }
00790
00791 bool DOFMapBuilder::regionIsMaximal(int r) const
00792 {
00793 const CellFilterStub* reg = fsr_->region(r).get();
00794 return (dynamic_cast<const MaximalCellFilter*>(reg) != 0);
00795 }
00796
00797 void DOFMapBuilder::markBCRows(int block)
00798 {
00799 isBCRow_[block] = rcp(new Array<int>(rowMap_[block]->numLocalDOFs()));
00800 int ndof = rowMap_[block]->numLocalDOFs();
00801 Array<int>& isBC = *isBCRow_[block];
00802 for (int i=0; i<ndof; i++) isBC[i] = false;
00803
00804 RCP<Array<int> > cellLID = rcp(new Array<int>());
00805 Array<RCP<Array<int> > > cellBatches;
00806 const RCP<DOFMapBase>& rowMap = rowMap_[block];
00807
00808 for (int r=0; r<fsr_->numRegions(); r++)
00809 {
00810
00811 CellFilter region = fsr_->region(r);
00812
00813 if (!fsr_->isBCRegion(r)) continue;
00814
00815 int dim = region.dimension(mesh_);
00816 CellSet cells = region.getCells(mesh_);
00817 cellLID->resize(0);
00818 for (CellIterator c=cells.begin(); c != cells.end(); c++)
00819 {
00820 cellLID->append(*c);
00821 }
00822 if (cellLID->size() == 0) continue;
00823
00824
00825 const Set<int>& allBcFuncs = fsr_->bcVarsOnRegion(r);
00826
00827 Set<int> bcFuncs;
00828 for (Set<int>::const_iterator
00829 i=allBcFuncs.begin(); i != allBcFuncs.end(); i++)
00830 {
00831 if (block == fsr_->blockForVarID(*i))
00832 {
00833 bcFuncs.put(fsr_->reducedVarID(*i));
00834 }
00835 }
00836 if (bcFuncs.size()==0) continue;
00837 Array<int> bcFuncID = bcFuncs.elements();
00838
00839 Array<Array<int> > dofs;
00840 Array<int> nNodes;
00841
00842 RCP<const MapStructure> s
00843 = rowMap->getDOFsForCellBatch(dim, *cellLID, bcFuncs, dofs, nNodes,0);
00844 int offset = rowMap->lowestLocalDOF();
00845 int high = offset + rowMap->numLocalDOFs();
00846
00847 for (int c=0; c<cellLID->size(); c++)
00848 {
00849 for (int b=0; b< s->numBasisChunks(); b++)
00850 {
00851 int nFuncs = s->numFuncs(b);
00852 for (int n=0; n<nNodes[b]; n++)
00853 {
00854 for (int f=0; f<bcFuncID.size(); f++)
00855 {
00856 int chunk = s->chunkForFuncID(bcFuncID[f]);
00857 if (chunk != b) continue;
00858 int funcOffset = s->indexForFuncID(bcFuncID[f]);
00859 int dof = dofs[b][(c*nFuncs + funcOffset)*nNodes[b]+n];
00860 if (dof < offset || dof >= high) continue;
00861 (*isBCRow_[block])[dof-offset]=true;
00862 }
00863 }
00864 }
00865 }
00866 }
00867 }
00868
00869
00870
00871 void DOFMapBuilder::markBCCols(int block)
00872 {
00873 isBCCol_[block] = rcp(new Array<int>(colMap_[block]->numLocalDOFs()));
00874 int ndof = colMap_[block]->numLocalDOFs();
00875 Array<int>& isBC = *isBCCol_[block];
00876 for (int i=0; i<ndof; i++) isBC[i] = false;
00877
00878 RCP<Array<int> > cellLID = rcp(new Array<int>());
00879 Array<RCP<Array<int> > > cellBatches;
00880 const RCP<DOFMapBase>& colMap = colMap_[block];
00881
00882 for (int r=0; r<fsr_->numRegions(); r++)
00883 {
00884
00885 CellFilter region = fsr_->region(r);
00886
00887 if (!fsr_->isBCRegion(r)) continue;
00888
00889 int dim = region.dimension(mesh_);
00890 CellSet cells = region.getCells(mesh_);
00891 cellLID->resize(0);
00892 for (CellIterator c=cells.begin(); c != cells.end(); c++)
00893 {
00894 cellLID->append(*c);
00895 }
00896 if (cellLID->size() == 0) continue;
00897
00898
00899 const Set<int>& allBcFuncs = fsr_->bcUnksOnRegion(r);
00900
00901 Set<int> bcFuncs;
00902 for (Set<int>::const_iterator
00903 i=allBcFuncs.begin(); i != allBcFuncs.end(); i++)
00904 {
00905 if (block == fsr_->blockForUnkID(*i))
00906 {
00907 bcFuncs.put(fsr_->reducedUnkID(*i));
00908 }
00909 }
00910 if (bcFuncs.size()==0) continue;
00911 Array<int> bcFuncID = bcFuncs.elements();
00912
00913 Array<Array<int> > dofs;
00914 Array<int> nNodes;
00915
00916 RCP<const MapStructure> s
00917 = colMap->getDOFsForCellBatch(dim, *cellLID, bcFuncs, dofs, nNodes,0);
00918 int offset = colMap->lowestLocalDOF();
00919 int high = offset + colMap->numLocalDOFs();
00920
00921 for (int c=0; c<cellLID->size(); c++)
00922 {
00923 for (int b=0; b< s->numBasisChunks(); b++)
00924 {
00925 int nFuncs = s->numFuncs(b);
00926 for (int n=0; n<nNodes[b]; n++)
00927 {
00928 for (int f=0; f<bcFuncID.size(); f++)
00929 {
00930 int chunk = s->chunkForFuncID(bcFuncID[f]);
00931 if (chunk != b) continue;
00932 int funcOffset = s->indexForFuncID(bcFuncID[f]);
00933 int dof = dofs[b][(c*nFuncs + funcOffset)*nNodes[b]+n];
00934 if (dof < offset || dof >= high)
00935 {
00936 remoteBCCols_[block]->insert(dof);
00937 }
00938 else
00939 {
00940 (*isBCCol_[block])[dof-offset]=true;
00941 }
00942 }
00943 }
00944 }
00945 }
00946 }
00947 }
00948
00949
00950
00951 namespace Sundance
00952 {
00953
00954 Array<Array<BasisFamily> > testBasisArray(const RCP<FunctionSupportResolver>& fsr)
00955 {
00956 Array<Array<BasisFamily> > rtn(fsr->numVarBlocks());
00957 for (int b=0; b<fsr->numVarBlocks(); b++)
00958 {
00959 for (int i=0; i<fsr->numVarIDs(b); i++)
00960 {
00961 rtn[b].append(BasisFamily::getBasis(fsr->varFuncData(b, i)));
00962 }
00963 }
00964 return rtn;
00965 }
00966
00967
00968 Array<Array<BasisFamily> > unkBasisArray(const RCP<FunctionSupportResolver>& fsr)
00969 {
00970 Array<Array<BasisFamily> > rtn(fsr->numUnkBlocks());
00971 for (int b=0; b<fsr->numUnkBlocks(); b++)
00972 {
00973 for (int i=0; i<fsr->numUnkIDs(b); i++)
00974 {
00975 rtn[b].append(BasisFamily::getBasis(fsr->unkFuncData(b, i)));
00976 }
00977 }
00978 return rtn;
00979 }
00980
00981
00982 }
