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 "SundanceSubmaximalNodalDOFMap.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
00046 SubmaximalNodalDOFMap
00047 ::SubmaximalNodalDOFMap(const Mesh& mesh,
00048 const CellFilter& cf,
00049 int nFuncs,
00050 int setupVerb)
00051 : DOFMapBase(mesh, setupVerb),
00052 dim_(0),
00053 nTotalFuncs_(nFuncs),
00054 domain_(cf),
00055 domains_(tuple(cf)),
00056 nodeLIDs_(),
00057 nodeDOFs_(),
00058 lidToPtrMap_(),
00059 mapStructure_()
00060 {
00061 Tabs tab0(0);
00062 SUNDANCE_MSG1(setupVerb, tab0 << "in SubmaximalNodalDOFMap ctor");
00063 Tabs tab1;
00064 SUNDANCE_MSG2(setupVerb, tab1 << "domain " << domain_);
00065 SUNDANCE_MSG2(setupVerb, tab1 << "N funcs " << nFuncs);
00066
00067 const MPIComm& comm = mesh.comm();
00068 int rank = comm.getRank();
00069 int nProc = comm.getNProc();
00070
00071 dim_ = cf.dimension(mesh);
00072 TEUCHOS_TEST_FOR_EXCEPT(dim_ != 0);
00073
00074 CellSet nodes = cf.getCells(mesh);
00075 int nc = nodes.numCells();
00076 nodeLIDs_.reserve(nc);
00077 nodeDOFs_.reserve(nc);
00078
00079 Array<Array<int> > remoteNodes(nProc);
00080
00081 int nextDOF = 0;
00082 int k=0;
00083 for (CellIterator c=nodes.begin(); c!=nodes.end(); c++, k++)
00084 {
00085 int nodeLID = *c;
00086 lidToPtrMap_.put(nodeLID, k);
00087 nodeLIDs_.append(nodeLID);
00088 int remoteOwner = rank;
00089 if (isRemote(0, nodeLID, remoteOwner))
00090 {
00091 int GID = mesh.mapLIDToGID(0, nodeLID);
00092 remoteNodes[remoteOwner].append(GID);
00093 for (int f=0; f<nFuncs; f++) nodeDOFs_.append(-1);
00094 }
00095 else
00096 {
00097 for (int f=0; f<nFuncs; f++) nodeDOFs_.append(nextDOF++);
00098 }
00099 }
00100
00101
00102 int localCount = nextDOF;
00103 computeOffsets(localCount);
00104
00105
00106 shareRemoteDOFs(remoteNodes);
00107
00108 BasisFamily basis = new Lagrange(1);
00109 mapStructure_ = rcp(new MapStructure(nTotalFuncs_, basis.ptr()));
00110 }
00111
00112 void SubmaximalNodalDOFMap::computeOffsets(int localCount)
00113 {
00114 Array<int> dofOffsets;
00115 int totalDOFCount;
00116 int myOffset = 0;
00117 int np = mesh().comm().getNProc();
00118 if (np > 1)
00119 {
00120 MPIContainerComm<int>::accumulate(localCount, dofOffsets, totalDOFCount,
00121 mesh().comm());
00122 myOffset = dofOffsets[mesh().comm().getRank()];
00123
00124 int nDofs = nodeDOFs_.size();
00125 for (int i=0; i<nDofs; i++)
00126 {
00127 if (nodeDOFs_[i] >= 0) nodeDOFs_[i] += myOffset;
00128 }
00129 }
00130 else
00131 {
00132 totalDOFCount = localCount;
00133 }
00134
00135 setLowestLocalDOF(myOffset);
00136 setNumLocalDOFs(localCount);
00137 setTotalNumDOFs(totalDOFCount);
00138 }
00139
00140
00141 void SubmaximalNodalDOFMap::shareRemoteDOFs(
00142 const Array<Array<int> >& outgoingNodeRequests)
00143 {
00144 int np = mesh().comm().getNProc();
00145 if (np==1) return;
00146 int rank = mesh().comm().getRank();
00147
00148 int nFuncs = nTotalFuncs_;
00149
00150 Array<Array<int> > incomingNodeRequests;
00151 Array<Array<int> > outgoingDOFs(np);
00152 Array<Array<int> > incomingDOFs;
00153
00154 SUNDANCE_MSG2(setupVerb(),
00155 "p=" << mesh().comm().getRank()
00156 << "synchronizing DOFs for cells of dimension 0");
00157 SUNDANCE_MSG2(setupVerb(),
00158 "p=" << mesh().comm().getRank()
00159 << " sending cell reqs d=0, GID="
00160 << outgoingNodeRequests);
00161
00162
00163 MPIContainerComm<int>::allToAll(outgoingNodeRequests,
00164 incomingNodeRequests,
00165 mesh().comm());
00166
00167
00168
00169 for (int p=0; p<np; p++)
00170 {
00171 if (p==rank) continue;
00172 const Array<int>& requestsFromProc = incomingNodeRequests[p];
00173 int nReq = requestsFromProc.size();
00174
00175 SUNDANCE_MSG4(setupVerb(), "p=" << mesh().comm().getRank()
00176 << " recv'd from proc=" << p
00177 << " reqs for DOFs for cells "
00178 << requestsFromProc);
00179
00180 outgoingDOFs[p].resize(nReq);
00181
00182 for (int c=0; c<nReq; c++)
00183 {
00184 int GID = requestsFromProc[c];
00185 SUNDANCE_MSG3(setupVerb(),
00186 "p=" << rank
00187 << " processing zero-cell with GID=" << GID);
00188 int LID = mesh().mapGIDToLID(0, GID);
00189 SUNDANCE_MSG3(setupVerb(),
00190 "p=" << rank
00191 << " LID=" << LID << " dofs="
00192 << nodeDOFs_[LID*nFuncs]);
00193 outgoingDOFs[p][c] = nodeDOFs_[LID*nFuncs];
00194 SUNDANCE_MSG3(setupVerb(),
00195 "p=" << rank
00196 << " done processing cell with GID=" << GID);
00197 }
00198 }
00199
00200 SUNDANCE_MSG2(setupVerb(),
00201 "p=" << mesh().comm().getRank()
00202 << "answering DOF requests for cells of dimension 0");
00203
00204
00205 MPIContainerComm<int>::allToAll(outgoingDOFs,
00206 incomingDOFs,
00207 mesh().comm());
00208
00209 SUNDANCE_MSG2(setupVerb(),
00210 "p=" << mesh().comm().getRank()
00211 << "communicated DOF answers for cells of dimension 0" );
00212
00213
00214
00215
00216 for (int p=0; p<mesh().comm().getNProc(); p++)
00217 {
00218 if (p==mesh().comm().getRank()) continue;
00219 const Array<int>& dofsFromProc = incomingDOFs[p];
00220 int numCells = dofsFromProc.size();
00221 for (int c=0; c<numCells; c++)
00222 {
00223 int cellGID = outgoingNodeRequests[p][c];
00224 int cellLID = mesh().mapGIDToLID(0, cellGID);
00225 int dof = dofsFromProc[c];
00226 for (int i=0; i<nFuncs; i++)
00227 {
00228 nodeDOFs_[cellLID*nFuncs + i] = dof+i;
00229 addGhostIndex(dof+i);
00230 }
00231 }
00232 }
00233 }
00234
00235
00236
00237 RCP<const Set<int> >
00238 SubmaximalNodalDOFMap::allowedFuncsOnCellBatch(int cellDim,
00239 const Array<int>& cellLID) const
00240 {
00241 Set<int> rtn;
00242 for (int f=0; f<nTotalFuncs_; f++) rtn.put(f);
00243 return rcp(new Set<int>(rtn));
00244 }
00245
00246
00247 RCP<const MapStructure>
00248 SubmaximalNodalDOFMap::getDOFsForCellBatch(int cellDim,
00249 const Array<int>& cellLID,
00250 const Set<int>& requestedFuncSet,
00251 Array<Array<int> >& dofs,
00252 Array<int>& nNodes,
00253 int verb) const
00254 {
00255 TimeMonitor timer(batchedDofLookupTimer());
00256 Tabs tab0;
00257
00258 SUNDANCE_MSG2(verb, tab0 << "in SubmaximalNodalDOFMap::getDOFsForCellBatch()");
00259
00260 TEUCHOS_TEST_FOR_EXCEPT(cellDim != 0);
00261
00262 dofs.resize(1);
00263 nNodes.resize(1);
00264 nNodes[0] = 1;
00265
00266 dofs[0].reserve(cellLID.size()*nTotalFuncs_);
00267
00268 for (int c=0; c<cellLID.size(); c++)
00269 {
00270 int lid = cellLID[c];
00271 int ptr = lidToPtrMap_.get(lid);
00272 for (int f=0; f<nTotalFuncs_; f++)
00273 {
00274 int q = ptr*nTotalFuncs_;
00275 dofs[0].append(nodeDOFs_[q + f]);
00276 }
00277 }
00278
00279 return mapStructure_;
00280 }
00281
00282
00283
00284 void SubmaximalNodalDOFMap::print(std::ostream& os) const
00285 {
00286 Tabs tab0(0);
00287 Out::os() << tab0 << "submaximal nodal dof map: " << std::endl;
00288 Tabs tab1;
00289 Out::os() << tab1 << "0-cells: " << std::endl;
00290 for (int c=0; c<nodeLIDs_.size(); c++)
00291 {
00292 Tabs tab2;
00293 int gid = mesh().mapLIDToGID(0, nodeLIDs_[c]);
00294 Out::os() << tab2 << "G=" << gid << " dofs={";
00295 for (int f=0; f<nTotalFuncs_; f++)
00296 {
00297 if (f != 0) Out::os() << ", ";
00298 Out::os() << nodeDOFs_[c*nTotalFuncs_ + f];
00299 }
00300 Out::os() << "}" << std::endl;
00301 }
00302 }
