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 "SundanceMap.hpp"
00032 #include "PlayaTabs.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "SundanceOrderedTuple.hpp"
00035 #include "SundancePartialElementDOFMap.hpp"
00036 #include "SundanceLagrange.hpp"
00037 #include "PlayaMPIContainerComm.hpp"
00038 #include "Teuchos_TimeMonitor.hpp"
00039
00040 using namespace Sundance;
00041 using namespace Teuchos;
00042 using Playa::MPIComm;
00043 using Playa::MPIContainerComm;
00044
00045 PartialElementDOFMap::PartialElementDOFMap(const Mesh& mesh,
00046 const CellFilter& subdomain,
00047 int nFuncs,
00048 int setupVerb)
00049 : DOFMapBase(mesh, setupVerb),
00050 dim_(mesh.spatialDim()),
00051 nFuncs_(nFuncs),
00052 nElems_(mesh.numCells(mesh.spatialDim())),
00053 subdomain_(subdomain),
00054 funcDomains_(nFuncs, subdomain),
00055 elemDofs_(nElems_ * nFuncs, -1),
00056 structure_(rcp(new MapStructure(nFuncs_, rcp(new Lagrange(0))))),
00057 allFuncs_()
00058 {
00059 init();
00060 Set<int> tmp;
00061 for (int i=0; i<nFuncs_; i++) tmp.put(i);
00062 allFuncs_ = rcp(new Set<int>(tmp));
00063 }
00064
00065
00066 RCP<const MapStructure>
00067 PartialElementDOFMap::getDOFsForCellBatch(int cellDim,
00068 const Array<int>& cellLID,
00069 const Set<int>& requestedFuncSet,
00070 Array<Array<int> >& dofs,
00071 Array<int>& nNodes,
00072 int verbosity) const
00073 {
00074 TimeMonitor timer(batchedDofLookupTimer());
00075
00076
00077 Tabs tab;
00078 SUNDANCE_MSG3(verbosity,
00079 tab << "PartialElementDOFMap::getDOFsForCellBatch(): cellDim="
00080 << cellDim
00081 << " cellLID=" << cellLID);
00082
00083
00084 dofs.resize(1);
00085 nNodes.resize(1);
00086 nNodes[0] = 1;
00087
00088 int nCells = cellLID.size();
00089
00090
00091 if (cellDim == dim_)
00092 {
00093 dofs[0].resize(nCells * nFuncs_);
00094 Array<int>& dof0 = dofs[0];
00095
00096 for (int c=0; c<nCells; c++)
00097 {
00098 for (int i=0; i<nFuncs_; i++)
00099 {
00100 dof0[c*nFuncs_ + i] = elemDofs_[cellLID[c]*nFuncs_+i];
00101 }
00102 }
00103 }
00104 else
00105 {
00106 dofs[0].resize(nCells * nFuncs_);
00107 Array<int> cofacetLIDs(nCells);
00108
00109 for (int c=0; c<nCells; c++)
00110 {
00111 TEUCHOS_TEST_FOR_EXCEPTION(mesh().numMaxCofacets(cellDim, cellLID[c]) > 1,
00112 std::runtime_error,
00113 "Attempt to do a trace of a L0 basis on a "
00114 "lower-dimensional cell having more than one "
00115 "maximal cofacets");
00116 int myFacetIndex = -1;
00117 cofacetLIDs[c] = mesh().maxCofacetLID(cellDim, cellLID[c],
00118 0, myFacetIndex);
00119 }
00120
00121 getDOFsForCellBatch(dim_, cofacetLIDs, requestedFuncSet, dofs,
00122 nNodes, verbosity);
00123 }
00124
00125 return structure_;
00126 }
00127
00128
00129
00130 void PartialElementDOFMap::init()
00131 {
00132 Tabs tab;
00133
00134 SUNDANCE_MSG1(setupVerb(), tab << "initializing element DOF map");
00135
00136 int cellDim = mesh().spatialDim();
00137
00138 Array<Array<int> > remoteElems(mesh().comm().getNProc());
00139
00140
00141 int nextDOF = 0;
00142 int owner;
00143
00144
00145 CellSet maxCells = subdomain_.getCells(mesh());
00146
00147 for (CellIterator iter=maxCells.begin(); iter != maxCells.end(); iter++)
00148 {
00149 int cellLID = *iter;
00150
00151 if (isRemote(cellDim, cellLID, owner))
00152 {
00153 int elemGID = mesh().mapLIDToGID(cellDim, cellLID);
00154 remoteElems[owner].append(elemGID);
00155 }
00156 else
00157 {
00158
00159 for (int i=0; i<nFuncs_; i++)
00160 {
00161 elemDofs_[cellLID*nFuncs_ + i] = nextDOF;
00162 nextDOF++;
00163 }
00164 }
00165 }
00166
00167
00168 int localCount = nextDOF;
00169 computeOffsets(localCount);
00170
00171
00172 shareRemoteDOFs(remoteElems);
00173 }
00174
00175
00176 RCP<const Set<int> > PartialElementDOFMap
00177 ::allowedFuncsOnCellBatch(int cellDim,
00178 const Array<int>& cellLID) const
00179 {
00180 static RCP<const Set<int> > empty = rcp(new Set<int>());
00181
00182 if (cellDim != dim_) return empty;
00183 for (int i=0; i<cellLID.size(); i++)
00184 {
00185 if (elemDofs_[cellLID[i]*nFuncs_] < 0) return empty;
00186 }
00187 return allFuncs_;
00188 }
00189
00190
00191 void PartialElementDOFMap::computeOffsets(int localCount)
00192 {
00193 Array<int> dofOffsets;
00194 int totalDOFCount;
00195 int myOffset = 0;
00196 int np = mesh().comm().getNProc();
00197 if (np > 1)
00198 {
00199 MPIContainerComm<int>::accumulate(localCount, dofOffsets, totalDOFCount,
00200 mesh().comm());
00201 myOffset = dofOffsets[mesh().comm().getRank()];
00202
00203 int nDofs = nElems_ * nFuncs_;
00204 for (int i=0; i<nDofs; i++)
00205 {
00206 if (elemDofs_[i] >= 0) elemDofs_[i] += myOffset;
00207 }
00208 }
00209 else
00210 {
00211 totalDOFCount = localCount;
00212 }
00213
00214 setLowestLocalDOF(myOffset);
00215 setNumLocalDOFs(localCount);
00216 setTotalNumDOFs(totalDOFCount);
00217 }
00218
00219
00220 void PartialElementDOFMap::shareRemoteDOFs(const Array<Array<int> >& outgoingCellRequests)
00221 {
00222 int cellDim = mesh().spatialDim();
00223
00224 int np = mesh().comm().getNProc();
00225 if (np==1) return;
00226 int rank = mesh().comm().getRank();
00227
00228 Array<Array<int> > incomingCellRequests;
00229 Array<Array<int> > outgoingDOFs(np);
00230 Array<Array<int> > incomingDOFs;
00231
00232 SUNDANCE_MSG2(setupVerb(),
00233 "p=" << mesh().comm().getRank()
00234 << "synchronizing DOFs for cells of dimension 0");
00235 SUNDANCE_MSG2(setupVerb(),
00236 "p=" << mesh().comm().getRank()
00237 << " sending cell reqs d=0, GID="
00238 << outgoingCellRequests);
00239
00240
00241 MPIContainerComm<int>::allToAll(outgoingCellRequests,
00242 incomingCellRequests,
00243 mesh().comm());
00244
00245
00246
00247 for (int p=0; p<np; p++)
00248 {
00249 if (p==rank) continue;
00250 const Array<int>& requestsFromProc = incomingCellRequests[p];
00251 int nReq = requestsFromProc.size();
00252
00253 SUNDANCE_MSG4(setupVerb(), "p=" << mesh().comm().getRank()
00254 << " recv'd from proc=" << p
00255 << " reqs for DOFs for cells "
00256 << requestsFromProc);
00257
00258 outgoingDOFs[p].resize(nReq);
00259
00260 for (int c=0; c<nReq; c++)
00261 {
00262 int GID = requestsFromProc[c];
00263 SUNDANCE_MSG3(setupVerb(),
00264 "p=" << rank
00265 << " processing zero-cell with GID=" << GID);
00266 int LID = mesh().mapGIDToLID(cellDim, GID);
00267 SUNDANCE_MSG3(setupVerb(),
00268 "p=" << rank
00269 << " LID=" << LID << " dofs="
00270 << elemDofs_[LID*nFuncs_]);
00271 outgoingDOFs[p][c] = elemDofs_[LID*nFuncs_];
00272 SUNDANCE_MSG3(setupVerb(),
00273 "p=" << rank
00274 << " done processing cell with GID=" << GID);
00275 }
00276 }
00277
00278 SUNDANCE_MSG2(setupVerb(),
00279 "p=" << mesh().comm().getRank()
00280 << "answering DOF requests for cells of dimension 0");
00281
00282
00283 MPIContainerComm<int>::allToAll(outgoingDOFs,
00284 incomingDOFs,
00285 mesh().comm());
00286
00287 SUNDANCE_MSG2(setupVerb(),
00288 "p=" << mesh().comm().getRank()
00289 << "communicated DOF answers for cells of dimension 0" );
00290
00291
00292
00293
00294 for (int p=0; p<mesh().comm().getNProc(); p++)
00295 {
00296 if (p==mesh().comm().getRank()) continue;
00297 const Array<int>& dofsFromProc = incomingDOFs[p];
00298 int numCells = dofsFromProc.size();
00299 for (int c=0; c<numCells; c++)
00300 {
00301 int cellGID = outgoingCellRequests[p][c];
00302 int cellLID = mesh().mapGIDToLID(cellDim, cellGID);
00303 int dof = dofsFromProc[c];
00304 for (int i=0; i<nFuncs_; i++)
00305 {
00306 elemDofs_[cellLID*nFuncs_ + i] = dof+i;
00307 addGhostIndex(dof+i);
00308 }
00309 }
00310 }
00311 }
00312
00313
00314 void PartialElementDOFMap::print(std::ostream& os) const
00315 {
00316 os << elemDofs_ << std::endl;
00317 }
