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 "SundanceInhomogeneousEdgeLocalizedDOFMap.hpp"
00032
00033 #include "SundanceEdgeLocalizedBasis.hpp"
00034 #include "SundanceOrderedTuple.hpp"
00035
00036 #include <numeric>
00037 #include <algorithm>
00038
00039 namespace Sundance {
00040
00041 using Teuchos::null;
00042
00043 InhomogeneousEdgeLocalizedDOFMap::InhomogeneousEdgeLocalizedDOFMap(const Mesh& mesh,
00044 const Array<Map<Set<int>, CellFilter> >& funcSetToDomainMap,
00045 int setupVerb) :
00046 DOFMapBase(mesh, setupVerb),
00047 funcDomains_(),
00048 edgeDofs_()
00049 {
00050 SUNDANCE_MSG1(setupVerb, "in InhomogeneousEdgeLocalizedDOFMap ctor");
00051
00052 TEUCHOS_TEST_FOR_EXCEPTION(mesh.comm().getNProc() != 1,
00053 std::runtime_error,
00054 "distributed inhomogeneous edge localized DOF maps not yet supported");
00055
00056 SUNDANCE_MSG4(setupVerb, "func set to domain map " << funcSetToDomainMap);
00057
00058 TEUCHOS_ASSERT(funcSetToDomainMap.size() > 1);
00059 TEUCHOS_ASSERT(funcSetToDomainMap[0].empty());
00060
00061 TEUCHOS_ASSERT(funcSetToDomainMap.size() <= this->meshDimension() + 1);
00062
00063 Map<int, Array<Array<CellFilter> > > funcFilters;
00064 Map<CellFilter, Array<int> > filterEdges;
00065 for (int dim = 1; dim <= this->meshDimension(); ++dim) {
00066 const Map<Set<int>, CellFilter> &funcsOnDom = funcSetToDomainMap[dim];
00067 for (Map<Set<int>, CellFilter>::const_iterator it = funcsOnDom.begin(),
00068 it_end = funcsOnDom.end();
00069 it != it_end;
00070 ++it)
00071 {
00072 const CellFilter &subdomain = it->second;
00073 filterEdges[subdomain] = getEdgeLIDs(subdomain);
00074
00075 const Set<int> &funcs = it->first;
00076 for (Set<int>::const_iterator jt = funcs.begin(),
00077 jt_end = funcs.end();
00078 jt != jt_end;
00079 ++jt)
00080 {
00081 Array<Array<CellFilter> > &filters = funcFilters[*jt];
00082 filters.resize(this->meshDimension() + 1);
00083 filters[dim].append(subdomain);
00084 }
00085 }
00086 }
00087
00088 SUNDANCE_MSG4(setupVerb, "subdomains to edges " << filterEdges);
00089 SUNDANCE_MSG4(setupVerb, "func ids to subdomains " << funcFilters);
00090
00091 const int functionCount = funcFilters.empty() ? 0 : funcFilters.rbegin()->first + 1;
00092 funcDomains_.resize(functionCount);
00093
00094 for (Map<int, Array<Array<CellFilter> > >::const_iterator funcIt = funcFilters.begin(),
00095 funcItEnd = funcFilters.end();
00096 funcIt != funcItEnd;
00097 ++funcIt)
00098 {
00099 const Array<CellFilter> &maxDimFilters = (funcIt->second)[this->meshDimension()];
00100 funcDomains_[funcIt->first] = std::accumulate(maxDimFilters.begin(), maxDimFilters.end(), CellFilter());
00101 }
00102
00103 SUNDANCE_MSG4(setupVerb, "func ids to max dim subdomains " << funcDomains_);
00104
00105 typedef OrderedPair<int, int> DofId;
00106 Array<DofId> dofSet;
00107 for (Map<int, Array<Array<CellFilter> > >::const_iterator funcIt = funcFilters.begin(),
00108 funcItEnd = funcFilters.end();
00109 funcIt !=funcItEnd;
00110 ++funcIt)
00111 {
00112 const int funcId = funcIt->first;
00113 Set<int> dofBearingEdges;
00114 for (int dim = 1; dim <= this->meshDimension(); ++dim)
00115 {
00116 const Array<CellFilter> &dimFilters = (funcIt->second)[dim];
00117 for (Array<CellFilter>::const_iterator filtIt = dimFilters.begin(),
00118 filtItEnd = dimFilters.end();
00119 filtIt != filtItEnd;
00120 ++filtIt)
00121 {
00122 const Array<int> &edges = filterEdges[*filtIt];
00123 dofBearingEdges.insert(edges.begin(), edges.end());
00124 }
00125 }
00126 for (Set<int>::const_iterator edgeIt = dofBearingEdges.begin(),
00127 edgeItEnd = dofBearingEdges.end();
00128 edgeIt != edgeItEnd;
00129 ++edgeIt)
00130 {
00131 dofSet.push_back(DofId(*edgeIt, funcId));
00132 }
00133 }
00134
00135 const int dofCount = dofSet.size();
00136
00137 SUNDANCE_MSG1(setupVerb, "number of degrees of freedom = " << dofCount);
00138 SUNDANCE_MSG4(setupVerb, "set of degrees of freedom " << dofSet);
00139
00140 const int edgeCount = mesh.numCells(1);
00141 edgeDofs_.resize(edgeCount, Array<int>(functionCount, -1));
00142
00143 std::sort(dofSet.begin(), dofSet.end());
00144 int dofRank = 0;
00145 for (Array<DofId>::const_iterator dofIt = dofSet.begin(),
00146 dofItEnd = dofSet.end();
00147 dofIt != dofItEnd;
00148 ++dofIt)
00149 {
00150 const int edgeId = dofIt->first();
00151 const int funcId = dofIt->second();
00152 edgeDofs_[edgeId][funcId] = dofRank++;
00153 }
00154
00155 SUNDANCE_MSG4(setupVerb, "degrees of freedom on edges " << edgeDofs_);
00156
00157 this->setLowestLocalDOF(0);
00158 this->setNumLocalDOFs(dofCount);
00159 this->setTotalNumDOFs(dofCount);
00160 }
00161
00162 RCP<const MapStructure>
00163 InhomogeneousEdgeLocalizedDOFMap::getDOFsForCellBatch(int cellDim,
00164 const Array<int>& cellLIDs,
00165 const Set<int>& requestedFuncSet,
00166 Array<Array<int> >& dofs,
00167 Array<int>& nNodes,
00168 int verb) const
00169 {
00170 TEUCHOS_ASSERT(requestedFuncSet.setDifference(*this->allowedFuncsOnCellBatch(cellDim, cellLIDs)).empty());
00171
00172 int edgesPerCell;
00173 if (cellDim == 0)
00174 {
00175
00176 edgesPerCell = 0;
00177 const Array<int> empty;
00178 getDOFsForEdgeBatch(empty, requestedFuncSet, dofs, verb);
00179 }
00180 else if (cellDim == 1)
00181 {
00182
00183 edgesPerCell = 1;
00184 getDOFsForEdgeBatch(cellLIDs, requestedFuncSet, dofs, verb);
00185 }
00186 else
00187 {
00188
00189 const int edgeDim = 1;
00190 edgesPerCell = this->mesh().numFacets(cellDim, 0, edgeDim);
00191 Array<int> edgeLIDs;
00192 {
00193 Array<int> dummyOrientations;
00194 this->mesh().getFacetLIDs(cellDim, cellLIDs, edgeDim, edgeLIDs, dummyOrientations);
00195 TEUCHOS_ASSERT(edgeLIDs.size() == edgesPerCell * cellLIDs.size());
00196 }
00197
00198 getDOFsForEdgeBatch(edgeLIDs, requestedFuncSet, dofs, verb);
00199 }
00200
00201 nNodes.resize(1);
00202 nNodes[0] = edgesPerCell;
00203
00204 const Array<int> requestedFuncArray(requestedFuncSet.begin(), requestedFuncSet.end());
00205 return rcp(new MapStructure(requestedFuncSet.size(), rcp(new EdgeLocalizedBasis()), tuple(requestedFuncArray)));
00206 }
00207
00208 void
00209 InhomogeneousEdgeLocalizedDOFMap::getDOFsForEdgeBatch(const Array<int>& edgeLIDs,
00210 const Set<int>& requestedFuncSet,
00211 Array<Array<int> >& dofs,
00212 int verb) const
00213 {
00214 dofs.resize(1);
00215 Array<int> &dofChunk = dofs[0];
00216 dofChunk.resize(edgeLIDs.size() * requestedFuncSet.size());
00217
00218 Array<int>::iterator entryIt = dofChunk.begin();
00219 for (Array<int>::const_iterator edgeIt = edgeLIDs.begin(),
00220 edgeItEnd = edgeLIDs.end();
00221 edgeIt != edgeItEnd;
00222 ++edgeIt)
00223 {
00224 for (Set<int>::const_iterator funcIt = requestedFuncSet.begin(),
00225 funcItEnd = requestedFuncSet.end();
00226 funcIt != funcItEnd;
00227 ++funcIt)
00228 {
00229 (*entryIt++) = edgeDofs_[*edgeIt][*funcIt];
00230 }
00231 }
00232 }
00233
00234 RCP<const Set<int> >
00235 InhomogeneousEdgeLocalizedDOFMap::allowedFuncsOnCellBatch(int cellDim,
00236 const Array<int>& cellLIDs) const
00237 {
00238 if (cellDim == 0)
00239 {
00240
00241 return rcp(new Set<int>());
00242 }
00243
00244 if (cellDim == this->meshDimension())
00245 {
00246
00247
00248 return this->allFuncIDs();
00249 }
00250
00251 const int edgeDim = 1;
00252 if (cellDim == edgeDim)
00253 {
00254
00255 return this->allowedFuncsOnEdgeBatch(cellLIDs);
00256 }
00257
00258
00259 Array<int> edgeLIDs;
00260 Array<int> dummyOrientations;
00261 this->mesh().getFacetLIDs(cellDim, cellLIDs, edgeDim, edgeLIDs, dummyOrientations);
00262 return this->allowedFuncsOnEdgeBatch(edgeLIDs);
00263 }
00264
00265 RCP<Set<int> >
00266 InhomogeneousEdgeLocalizedDOFMap::allowedFuncsOnEdgeBatch(const Array<int> &edgeLIDs) const
00267 {
00268 const RCP<Set<int> > result = this->allFuncIDs();
00269
00270 for (Array<int>::const_iterator edgeIt = edgeLIDs.begin(),
00271 edgeItEnd = edgeLIDs.end();
00272 edgeIt != edgeItEnd;
00273 ++edgeIt)
00274 {
00275 const Array<int> &dofs = edgeDofs_[*edgeIt];
00276
00277 for (Set<int>::iterator it = result->begin(); it != result->end(); )
00278 {
00279 if (dofs[*it] < 0)
00280 {
00281 result->erase(it++);
00282 }
00283 else
00284 {
00285 ++it;
00286 }
00287 }
00288
00289 if (result->empty())
00290 {
00291 break;
00292 }
00293 }
00294
00295 return result;
00296 }
00297
00298 RCP<Set<int> >
00299 InhomogeneousEdgeLocalizedDOFMap::allFuncIDs() const
00300 {
00301 const RCP<Set<int> > result(new Set<int>());
00302 for (int i = 0; i < funcDomains_.size(); ++i)
00303 {
00304 result->insert(result->end(), i);
00305 }
00306 return result;
00307 }
00308
00309 void
00310 InhomogeneousEdgeLocalizedDOFMap::print(std::ostream& os) const
00311 {
00312 os << "InhomogeneousEdgeLocalizedDOFMap\n";
00313 os << "Edges in mesh = " << edgeDofs_.size() << "\n";
00314 os << "Functions = " << funcDomains_.size() << "\n";
00315
00316 for (int dim = 1; dim <= this->meshDimension(); ++dim)
00317 {
00318 os << "Dofs on cells of dimension " << dim << ":\n";
00319 const int cellCount = this->mesh().numCells(dim);
00320 for (int iCell = 0; iCell < cellCount; ++iCell)
00321 {
00322 os << " Cell LID " << iCell << "\n";
00323 const RCP<const Set<int> > funcs = this->allowedFuncsOnCellBatch(dim, tuple(iCell));
00324 for (Set<int>::const_iterator it = funcs->begin(),
00325 it_end = funcs->end();
00326 it != it_end;
00327 ++it)
00328 {
00329 const int funcId = *it;
00330 os << " Function " << funcId << " : {";
00331 Array<int> dofs;
00332 this->getDOFsForCell(dim, iCell, funcId, dofs);
00333 for (Array<int>::const_iterator jt = dofs.begin(),
00334 jt_end = dofs.end();
00335 jt != jt_end;
00336 ++jt)
00337 {
00338 if (jt != dofs.begin())
00339 {
00340 os << ", ";
00341 }
00342 os << *jt;
00343 }
00344 os << "}\n";
00345 }
00346 }
00347 }
00348 }
00349
00350 int
00351 InhomogeneousEdgeLocalizedDOFMap::meshDimension() const
00352 {
00353 return this->mesh().spatialDim();
00354 }
00355
00356 Array<int>
00357 InhomogeneousEdgeLocalizedDOFMap::getEdgeLIDs(const CellFilter &filter) const
00358 {
00359 const int cellDim = filter.dimension(this->mesh());
00360
00361 const int nodeDim = 0;
00362 if (cellDim == nodeDim)
00363 {
00364
00365 return Array<int>();
00366 }
00367
00368 const CellSet cells = filter.getCells(this->mesh());
00369 const Array<int> cellLIDs(cells.begin(), cells.end());
00370
00371 const int edgeDim = 1;
00372 if (cellDim == edgeDim)
00373 {
00374
00375 return cellLIDs;
00376 }
00377
00378
00379 Array<int> edgeLIDs;
00380 Array<int> dummyOrientations;
00381 this->mesh().getFacetLIDs(cellDim, cellLIDs, edgeDim, edgeLIDs, dummyOrientations);
00382
00383 return edgeLIDs;
00384 }
00385
00386 }