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 "SundanceNodalDOFMap.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 NodalDOFMap::NodalDOFMap(const Mesh& mesh,
00046 int nFuncs,
00047 const CellFilter& maxCellFilter,
00048 int setupVerb)
00049 : SpatiallyHomogeneousDOFMapBase(mesh, nFuncs, setupVerb),
00050 maxCellFilter_(maxCellFilter),
00051 dim_(mesh.spatialDim()),
00052 nFuncs_(nFuncs),
00053 nElems_(mesh.numCells(mesh.spatialDim())),
00054 nNodes_(mesh.numCells(0)),
00055 nNodesPerElem_(mesh.numFacets(mesh.spatialDim(),0,0)),
00056 elemDofs_(nElems_ * nFuncs * nNodesPerElem_),
00057 nodeDofs_(mesh.numCells(0)*nFuncs_, -1),
00058 structure_(rcp(new MapStructure(nFuncs_, rcp(new Lagrange(1)))))
00059 {
00060 init();
00061 }
00062
00063
00064 RCP<const MapStructure>
00065 NodalDOFMap::getDOFsForCellBatch(int cellDim,
00066 const Array<int>& cellLID,
00067 const Set<int>& requestedFuncSet,
00068 Array<Array<int> >& dofs,
00069 Array<int>& nNodes,
00070 int verbosity) const
00071 {
00072 TimeMonitor timer(batchedDofLookupTimer());
00073
00074 Tabs tab;
00075 SUNDANCE_MSG3(verbosity,
00076 tab << "NodalDOFMap::getDOFsForCellBatch(): cellDim=" << cellDim
00077 << " cellLID=" << cellLID);
00078
00079
00080 dofs.resize(1);
00081 nNodes.resize(1);
00082
00083 int nCells = cellLID.size();
00084
00085
00086 if (cellDim == dim_)
00087 {
00088 int dofsPerElem = nFuncs_ * nNodesPerElem_;
00089 nNodes[0] = nNodesPerElem_;
00090 dofs[0].resize(nCells * dofsPerElem);
00091 Array<int>& dof0 = dofs[0];
00092
00093 for (int c=0; c<nCells; c++)
00094 {
00095 for (int i=0; i<dofsPerElem; i++)
00096 {
00097 dof0[c*dofsPerElem + i] = elemDofs_[cellLID[c]*dofsPerElem+i];
00098 }
00099 }
00100 }
00101 else if (cellDim == 0)
00102 {
00103 nNodes[0] = 1;
00104 dofs[0].resize(nCells * nFuncs_);
00105 Array<int>& dof0 = dofs[0];
00106
00107 for (int c=0; c<nCells; c++)
00108 {
00109 for (int i=0; i<nFuncs_; i++)
00110 {
00111 dof0[c*nFuncs_ + i] = nodeDofs_[cellLID[c]*nFuncs_+i];
00112 }
00113 }
00114 }
00115 else
00116 {
00117 int nFacets = mesh().numFacets(cellDim, cellLID[0], 0);
00118 nNodes[0] = nFacets;
00119 int dofsPerCell = nFuncs_ * nNodes[0];
00120 dofs[0].resize(nCells * dofsPerCell);
00121 Array<int>& dof0 = dofs[0];
00122 Array<int> facetLIDs(nCells * nNodes[0]);
00123 Array<int> dummy(nCells * nNodes[0]);
00124 mesh().getFacetLIDs(cellDim, cellLID, 0, facetLIDs, dummy);
00125
00126 for (int c=0; c<nCells; c++)
00127 {
00128 for (int f=0; f<nFacets; f++)
00129 {
00130 int facetCellLID = facetLIDs[c*nFacets+f];
00131 for (int i=0; i<nFuncs_; i++)
00132 {
00133 dof0[(c*nFuncs_+i)*nFacets + f]
00134 = nodeDofs_[facetCellLID*nFuncs_+i];
00135 }
00136 }
00137 }
00138 }
00139
00140 return structure_;
00141 }
00142
00143
00144
00145 void NodalDOFMap::init()
00146 {
00147 Tabs tab;
00148
00149 SUNDANCE_MSG1(setupVerb(), tab << "initializing nodal DOF map");
00150
00151 Array<Array<int> > remoteNodes(mesh().comm().getNProc());
00152 Array<int> hasProcessedCell(nNodes_, 0);
00153
00154
00155 int nextDOF = 0;
00156 int owner;
00157
00158 int nFacets = mesh().numFacets(dim_, 0, 0);
00159 Array<int> elemLID(nElems_);
00160 Array<int> facetLID(nFacets*nElems_);
00161 Array<int> facetOrientations(nFacets*nElems_);
00162
00163
00164 CellSet maxCells = maxCellFilter_.getCells(mesh());
00165
00166 int cc = 0;
00167 for (CellIterator iter=maxCells.begin(); iter != maxCells.end(); iter++, cc++)
00168 {
00169 int c = *iter;
00170 elemLID[cc] = c;
00171 }
00172
00173 mesh().getFacetLIDs(dim_, elemLID, 0, facetLID, facetOrientations);
00174
00175 for (int c=0; c<nElems_; c++)
00176 {
00177
00178 for (int f=0; f<nFacets; f++)
00179 {
00180
00181
00182 int fLID = facetLID[c*nFacets+f];
00183 if (hasProcessedCell[fLID] == 0)
00184 {
00185
00186 if (isRemote(0, fLID, owner))
00187 {
00188 int facetGID
00189 = mesh().mapLIDToGID(0, fLID);
00190 remoteNodes[owner].append(facetGID);
00191
00192 }
00193 else
00194 {
00195
00196 for (int i=0; i<nFuncs_; i++)
00197 {
00198 nodeDofs_[fLID*nFuncs_ + i] = nextDOF;
00199 nextDOF++;
00200 }
00201 }
00202 hasProcessedCell[fLID] = 1;
00203 }
00204 }
00205 }
00206
00207
00208
00209 int localCount = nextDOF;
00210 computeOffsets(localCount);
00211
00212
00213 shareRemoteDOFs(remoteNodes);
00214
00215
00216
00217 for (int c=0; c<nElems_; c++)
00218 {
00219
00220 for (int f=0; f<nFacets; f++)
00221 {
00222 int fLID = facetLID[c*nFacets+f];
00223 for (int i=0; i<nFuncs_; i++)
00224 {
00225 elemDofs_[(c*nFuncs_+i)*nFacets + f] = nodeDofs_[fLID*nFuncs_ + i];
00226 }
00227 }
00228 }
00229
00230 }
00231
00232
00233 void NodalDOFMap::computeOffsets(int localCount)
00234 {
00235 Array<int> dofOffsets;
00236 int totalDOFCount;
00237 int myOffset = 0;
00238 int np = mesh().comm().getNProc();
00239 if (np > 1)
00240 {
00241 MPIContainerComm<int>::accumulate(localCount, dofOffsets, totalDOFCount,
00242 mesh().comm());
00243 myOffset = dofOffsets[mesh().comm().getRank()];
00244
00245 int nDofs = nNodes_ * nFuncs_;
00246 for (int i=0; i<nDofs; i++)
00247 {
00248 if (nodeDofs_[i] >= 0) nodeDofs_[i] += myOffset;
00249 }
00250 }
00251 else
00252 {
00253 totalDOFCount = localCount;
00254 }
00255
00256 setLowestLocalDOF(myOffset);
00257 setNumLocalDOFs(localCount);
00258 setTotalNumDOFs(totalDOFCount);
00259 }
00260
00261
00262 void NodalDOFMap::shareRemoteDOFs(const Array<Array<int> >& outgoingCellRequests)
00263 {
00264 int np = mesh().comm().getNProc();
00265 if (np==1) return;
00266 int rank = mesh().comm().getRank();
00267
00268 Array<Array<int> > incomingCellRequests;
00269 Array<Array<int> > outgoingDOFs(np);
00270 Array<Array<int> > incomingDOFs;
00271
00272 SUNDANCE_MSG2(setupVerb(),
00273 "p=" << mesh().comm().getRank()
00274 << "synchronizing DOFs for cells of dimension 0");
00275 SUNDANCE_MSG2(setupVerb(),
00276 "p=" << mesh().comm().getRank()
00277 << " sending cell reqs d=0, GID="
00278 << outgoingCellRequests);
00279
00280
00281 MPIContainerComm<int>::allToAll(outgoingCellRequests,
00282 incomingCellRequests,
00283 mesh().comm());
00284
00285
00286
00287 for (int p=0; p<np; p++)
00288 {
00289 if (p==rank) continue;
00290 const Array<int>& requestsFromProc = incomingCellRequests[p];
00291 int nReq = requestsFromProc.size();
00292
00293 SUNDANCE_MSG4(setupVerb(), "p=" << mesh().comm().getRank()
00294 << " recv'd from proc=" << p
00295 << " reqs for DOFs for cells "
00296 << requestsFromProc);
00297
00298 outgoingDOFs[p].resize(nReq);
00299
00300 for (int c=0; c<nReq; c++)
00301 {
00302 int GID = requestsFromProc[c];
00303 SUNDANCE_MSG3(setupVerb(),
00304 "p=" << rank
00305 << " processing zero-cell with GID=" << GID);
00306 int LID = mesh().mapGIDToLID(0, GID);
00307 SUNDANCE_MSG3(setupVerb(),
00308 "p=" << rank
00309 << " LID=" << LID << " dofs="
00310 << nodeDofs_[LID*nFuncs_]);
00311 outgoingDOFs[p][c] = nodeDofs_[LID*nFuncs_];
00312 SUNDANCE_MSG3(setupVerb(),
00313 "p=" << rank
00314 << " done processing cell with GID=" << GID);
00315 }
00316 }
00317
00318 SUNDANCE_MSG2(setupVerb(),
00319 "p=" << mesh().comm().getRank()
00320 << "answering DOF requests for cells of dimension 0");
00321
00322
00323 MPIContainerComm<int>::allToAll(outgoingDOFs,
00324 incomingDOFs,
00325 mesh().comm());
00326
00327 SUNDANCE_MSG2(setupVerb(),
00328 "p=" << mesh().comm().getRank()
00329 << "communicated DOF answers for cells of dimension 0" );
00330
00331
00332
00333
00334 for (int p=0; p<mesh().comm().getNProc(); p++)
00335 {
00336 if (p==mesh().comm().getRank()) continue;
00337 const Array<int>& dofsFromProc = incomingDOFs[p];
00338 int numCells = dofsFromProc.size();
00339 for (int c=0; c<numCells; c++)
00340 {
00341 int cellGID = outgoingCellRequests[p][c];
00342 int cellLID = mesh().mapGIDToLID(0, cellGID);
00343 int dof = dofsFromProc[c];
00344 for (int i=0; i<nFuncs_; i++)
00345 {
00346 nodeDofs_[cellLID*nFuncs_ + i] = dof+i;
00347 addGhostIndex(dof+i);
00348 }
00349 }
00350 }
00351 }