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 "SundanceNodalDOFMapHN.hpp"
00036 #include "SundanceLagrange.hpp"
00037 #include "PlayaMPIContainerComm.hpp"
00038 #include "Teuchos_TimeMonitor.hpp"
00039
00040
00041 using namespace Sundance;
00042 using namespace Teuchos;
00043 using Playa::MPIComm;
00044 using Playa::MPIContainerComm;
00045
00046 NodalDOFMapHN::NodalDOFMapHN(const Mesh& mesh,
00047 int nFuncs,
00048 const CellFilter& maxCellFilter,
00049 int setupVerb)
00050 : HNDoFMapBaseHomogeneous(mesh, nFuncs, setupVerb),
00051 maxCellFilter_(maxCellFilter),
00052 dim_(mesh.spatialDim()),
00053 nFuncs_(nFuncs),
00054 nElems_(mesh.numCells(mesh.spatialDim())),
00055 nNodes_(mesh.numCells(0)),
00056 nNodesPerElem_(mesh.numFacets(mesh.spatialDim(),0,0)),
00057 elemDofs_(nElems_ * nFuncs * nNodesPerElem_),
00058 nodeDofs_(mesh.numCells(0)*nFuncs_, -2),
00059 structure_(rcp(new MapStructure(nFuncs_, rcp(new Lagrange(1))))),
00060 hasCellHanging_(nElems_),
00061 nodeIsHanging_(nElems_ * nFuncs * nNodesPerElem_),
00062 cellsWithHangingDoF_globalDoFs_(),
00063 cells_To_NodeLIDs_(),
00064 hangingNodeLID_to_NodesLIDs_(),
00065 hangindNodeLID_to_Coefs_(),
00066 matrixStore_(),
00067 facetLID_()
00068 {
00069
00070 matrixStore_.init(1);
00071 init();
00072 }
00073
00074 void NodalDOFMapHN::init()
00075 {
00076 Tabs tab;
00077
00078 SUNDANCE_MSG2(setupVerb(), tab << "NodalDOFMapHN initializing nodal DOF map nrFunc:"
00079 << nFuncs_ << " nNodes_:" << nNodes_ << " nElems_:" <<nElems_);
00080
00081 Array<Array<int> > remoteNodes(mesh().comm().getNProc());
00082 Array<int> hasProcessedCell(nNodes_, 0);
00083
00084
00085 int nextDOF = 0;
00086 int owner;
00087
00088 nFacets_ = mesh().numFacets(dim_, 0, 0);
00089 Array<int> elemLID(nElems_);
00090 Array<int> facetOrientations(nFacets_*nElems_);
00091
00092
00093 CellSet maxCells = maxCellFilter_.getCells(mesh());
00094
00095 int cc = 0;
00096 for (CellIterator iter=maxCells.begin(); iter != maxCells.end(); iter++, cc++)
00097 {
00098 int c = *iter;
00099 elemLID[cc] = c;
00100 }
00101
00102
00103 mesh().getFacetLIDs(dim_, elemLID, 0, facetLID_, facetOrientations);
00104
00105 for (int c=0; c<nElems_; c++) {
00106 hasCellHanging_[c] = false;
00107
00108 for (int f=0; f<nFacets_; f++) {
00109
00110
00111 int fLID = facetLID_[c*nFacets_+f];
00112 SUNDANCE_MSG2(setupVerb(), "NodalDOFMapHN::init() CellLID:"<< c <<"Try point LID:" << fLID << " facet:" << f);
00113 if (hasProcessedCell[fLID] == 0) {
00114
00115
00116 if ( isRemote(0, fLID, owner) && (!mesh().isElementHangingNode(0,fLID)) ) {
00117 int facetGID
00118 = mesh().mapLIDToGID(0, fLID);
00119 remoteNodes[owner].append(facetGID);
00120 }
00121 else {
00122 SUNDANCE_MSG2(setupVerb(), "NodalDOFMapHN::init() Doing point LID:" << fLID << " facet:" << f);
00123
00124 if ( mesh().isElementHangingNode(0,fLID) == false ){
00125
00126 for (int i=0; i<nFuncs_; i++){
00127 nodeDofs_[fLID*nFuncs_ + i] = nextDOF;
00128 nextDOF++;
00129 }
00130 } else {
00131 SUNDANCE_MSG2(setupVerb(), "NodalDOFMapHN::init() Hanging node found LID:" << fLID);
00132 hasCellHanging_[c] = true;
00133 for (int i=0; i<nFuncs_; i++){
00134 nodeDofs_[fLID*nFuncs_ + i] = -1;
00135 }
00136 Array<int> pointsLIDs;
00137 Array<int> facetIndex;
00138 Array<double> coefs;
00139
00140 getPointLIDsForHN( fLID, f, c , pointsLIDs, coefs , facetIndex);
00141
00142
00143 hangingNodeLID_to_NodesLIDs_.put( fLID , pointsLIDs );
00144 hangindNodeLID_to_Coefs_.put(fLID,coefs);
00145 }
00146 }
00147 hasProcessedCell[fLID] = 1;
00148 }
00149
00150 if (mesh().isElementHangingNode(0,fLID) == true) {
00151 hasCellHanging_[c] = true;
00152 }
00153 }
00154 }
00155
00156
00157
00158
00159 int localCount = nextDOF;
00160 computeOffsets(localCount);
00161
00162
00163 shareRemoteDOFs(remoteNodes);
00164
00165
00166
00167
00168 for (int c=0; c<nElems_; c++)
00169 {
00170 if (hasCellHanging_[c]){
00171 Array<int> HNDoFs(nFacets_);
00172 Array<double> transMatrix;
00173 Array<double> coefs;
00174
00175 transMatrix.resize(nFacets_*nFacets_);
00176 for (int ii = 0 ; ii < nFacets_*nFacets_ ; ii++) transMatrix[ii] = 0.0;
00177
00178 for (int f=0; f<nFacets_; f++)
00179 {
00180 int fLID = facetLID_[c*nFacets_+f];
00181 SUNDANCE_MSG2(setupVerb(), tab << "NodalDOFMapHN cell:" << c << " facetLID:" << fLID
00182 << " hanging:"<< nodeDofs_[fLID*nFuncs_] << " array:" << HNDoFs);
00183 if (nodeDofs_[fLID*nFuncs_] < 0)
00184 {
00185 Array<int> pointsLIDs;
00186 Array<int> facetIndex;
00187 Array<double> coefs;
00188
00189 getPointLIDsForHN( fLID, f, c , pointsLIDs, coefs , facetIndex);
00190
00191 for (int j=0 ; j < pointsLIDs.size() ; j++){
00192
00193 HNDoFs[facetIndex[j]] = pointsLIDs[j];
00194 transMatrix[f*nFacets_ + facetIndex[j]] = coefs[j];
00195 }
00196 }
00197 else
00198 {
00199
00200 HNDoFs[f] = fLID;
00201 transMatrix[f*nFacets_ + f] = 1.0;
00202 }
00203 }
00204
00205 int matrixIndex = matrixStore_.addMatrix(0,transMatrix);
00206 maxCellLIDwithHN_to_TrafoMatrix_.put( c , matrixIndex );
00207
00208
00209 cellsWithHangingDoF_globalDoFs_.put( c , HNDoFs );
00210
00211 SUNDANCE_MSG2(setupVerb(), tab << "NodalDOFMapHN initializing cellLID:" << c << " array:" << HNDoFs);
00212 SUNDANCE_MSG2(setupVerb(), tab << "NodalDOFMapHN initializing cellLID:" << c << " Trafo array:" << transMatrix);
00213
00214
00215 for (int f=0; f<nFacets_; f++)
00216 {
00217 int fLID = HNDoFs[f];
00218 for (int i=0; i<nFuncs_; i++) {
00219 elemDofs_[(c*nFuncs_+i)*nFacets_ + f] = nodeDofs_[fLID*nFuncs_ + i];
00220 }
00221
00222 }
00223 SUNDANCE_MSG2(setupVerb(),tab << "NodalDOFMapHN HN initializing cellLID:" << c << " elemDofs_:" << elemDofs_);
00224 }
00225 else {
00226
00227 for (int f=0; f<nFacets_; f++)
00228 {
00229 int fLID = facetLID_[c*nFacets_+f];
00230 for (int i=0; i<nFuncs_; i++)
00231 {
00232 elemDofs_[(c*nFuncs_+i)*nFacets_ + f] = nodeDofs_[fLID*nFuncs_ + i];
00233 }
00234 }
00235 SUNDANCE_MSG2(setupVerb(),tab << "NodalDOFMapHN initializing cellLID:" << c << " elemDofs_:" << elemDofs_);
00236 }
00237 }
00238 SUNDANCE_MSG2(setupVerb(),tab << "NodalDOFMapHN initializing DONE");
00239
00240 }
00241
00242 RCP<const MapStructure>
00243 NodalDOFMapHN::getDOFsForCellBatch(int cellDim,
00244 const Array<int>& cellLID,
00245 const Set<int>& requestedFuncSet,
00246 Array<Array<int> >& dofs,
00247 Array<int>& nNodes,
00248 int verbosity) const
00249 {
00250 TimeMonitor timer(batchedDofLookupTimer());
00251
00252 Tabs tab;
00253 SUNDANCE_MSG2(verbosity,
00254 tab << "NodalDOFMapHN::getDOFsForCellBatch(): cellDim=" << cellDim
00255 << " requestedFuncSet:" << requestedFuncSet
00256 << " cellLID=" << cellLID);
00257
00258
00259 dofs.resize(1);
00260 nNodes.resize(1);
00261
00262 int nCells = cellLID.size();
00263
00264
00265 if (cellDim == dim_)
00266 {
00267 int dofsPerElem = nFuncs_ * nNodesPerElem_;
00268 nNodes[0] = nNodesPerElem_;
00269 dofs[0].resize(nCells * dofsPerElem);
00270 Array<int>& dof0 = dofs[0];
00271 Array<int> tmpArray;
00272
00273 tmpArray.resize(dofsPerElem);
00274
00275 for (int c=0; c<nCells; c++)
00276 {
00277
00278 for (int i=0; i<dofsPerElem; i++) {
00279 dof0[c*dofsPerElem + i] = elemDofs_[cellLID[c]*dofsPerElem+i];
00280 tmpArray[i] = elemDofs_[cellLID[c]*dofsPerElem+i];
00281 }
00282 SUNDANCE_MSG2(verbosity,tab << "NodalDOFMapHN::getDOFsForCellBatch cellDim:" <<
00283 cellDim << " cellLID:" << c << " array:" << tmpArray);
00284 }
00285 }
00286 else if (cellDim == 0)
00287 {
00288 nNodes[0] = 1;
00289 dofs[0].resize(nCells * nFuncs_);
00290 Array<int>& dof0 = dofs[0];
00291 Array<int> tmpArray;
00292 tmpArray.resize(nFuncs_);
00293
00294 for (int c=0; c<nCells; c++)
00295 {
00296 for (int i=0; i<nFuncs_; i++)
00297 {
00298 dof0[c*nFuncs_ + i] = nodeDofs_[cellLID[c]*nFuncs_+i];
00299 tmpArray[i] = nodeDofs_[cellLID[c]*nFuncs_+i];
00300 }
00301 SUNDANCE_MSG2(verbosity,tab << "NodalDOFMapHN::getDOFsForCellBatch cellDim:" <<
00302 cellDim << " pointLID:" << cellLID[c] << " array:" << tmpArray);
00303 }
00304 }
00305 else
00306 {
00307
00308 int nFacets = mesh().numFacets(cellDim, cellLID[0], 0);
00309 nNodes[0] = nFacets;
00310 int dofsPerCell = nFuncs_ * nNodes[0];
00311 dofs[0].resize(nCells * dofsPerCell);
00312 Array<int>& dof0 = dofs[0];
00313 Array<int> facetLIDs(nCells * nNodes[0]);
00314 Array<int> dummy(nCells * nNodes[0]);
00315 Array<int> tmpArray;
00316 tmpArray.resize(nFuncs_*nFacets);
00317
00318 mesh().getFacetLIDs(cellDim, cellLID, 0, facetLIDs, dummy);
00319
00320 for (int c=0; c<nCells; c++)
00321 {
00322 for (int f=0; f<nFacets; f++)
00323 {
00324 int facetCellLID = facetLIDs[c*nFacets+f];
00325 for (int i=0; i<nFuncs_; i++)
00326 {
00327 dof0[(c*nFuncs_+i)*nFacets + f]
00328 = nodeDofs_[facetCellLID*nFuncs_+i];
00329
00330 if ( nodeDofs_[facetCellLID*nFuncs_+i] < 0){
00331 int tmp = 1;
00332
00333 int maxCell = mesh().maxCofacetLID( cellDim , cellLID[c] , 0 , tmp);
00334 Array<int> facetsLIDs;
00335 mesh().returnParentFacets( maxCell , 0 , facetsLIDs , tmp );
00336
00337 Array<int> childFacets(mesh().numFacets(mesh().spatialDim(),0,0));
00338 for (int jk = 0 ; jk < childFacets.size() ; jk++)
00339 childFacets[jk] = mesh().facetLID( mesh().spatialDim() , maxCell, 0 , jk , tmp);
00340 int fIndex = 0;
00341
00342 for (int jk = 0 ; jk < childFacets.size() ; jk++)
00343 if ( childFacets[jk] == facetCellLID) { fIndex = jk; break;}
00344 dof0[(c*nFuncs_+i)*nFacets + f] = nodeDofs_[facetsLIDs[fIndex]*nFuncs_ + i];
00345 }
00346 tmpArray[f*nFuncs_+i] = dof0[(c*nFuncs_+i)*nFacets + f];
00347 }
00348 }
00349 SUNDANCE_MSG2(verbosity,tab << "NodalDOFMapHN::getDOFsForCellBatch cellDim:" <<
00350 cellDim << " edge or face LID:" << cellLID[c] << " array:" << tmpArray);
00351 }
00352 }
00353 SUNDANCE_MSG2(verbosity,
00354 tab << "NodalDOFMapHN::getDOFsForCellBatch(): DONE");
00355 return structure_;
00356 }
00357
00358
00359 void NodalDOFMapHN::getTrafoMatrixForCell(
00360 int cellLID,
00361 int funcID,
00362 int& trafoMatrixSize,
00363 bool& doTransform,
00364 Array<double>& transfMatrix ) const {
00365
00366 trafoMatrixSize = nFacets_;
00367
00368 if (cellsWithHangingDoF_globalDoFs_.containsKey(cellLID))
00369 {
00370
00371 int matrixIndex = maxCellLIDwithHN_to_TrafoMatrix_.get( cellLID );
00372 matrixStore_.getMatrix( 0 , matrixIndex , transfMatrix );
00373 doTransform = true;
00374 }
00375 else
00376 {
00377 doTransform = false;
00378 }
00379 }
00380
00381 void NodalDOFMapHN::getTrafoMatrixForFacet(
00382 int cellDim,
00383 int cellLID,
00384 int facetIndex,
00385 int funcID,
00386 int& trafoMatrixSize,
00387 bool& doTransform,
00388 Array<double>& transfMatrix ) const {
00389
00390 int fIndex;
00391 int maxCellLID;
00392
00393 SUNDANCE_MSG2(setupVerb() , "NodalDOFMapHN::getTrafoMatrixForFacet() cellDim :" << cellDim << ", cellLID:" << cellLID);
00394 maxCellLID = mesh().maxCofacetLID( cellDim, cellLID, 0 , fIndex);
00395 SUNDANCE_MSG2(setupVerb() , "NodalDOFMapHN::getTrafoMatrixForFacet() testing :" << maxCellLID);
00396
00397
00398
00399 if (cellsWithHangingDoF_globalDoFs_.containsKey(maxCellLID))
00400 {
00401 int matrixIndex = maxCellLIDwithHN_to_TrafoMatrix_.get( maxCellLID );
00402 matrixStore_.getMatrix( 0 , matrixIndex , transfMatrix );
00403 doTransform = true;
00404 }
00405 else
00406 {
00407 doTransform = false;
00408 }
00409 }
00410
00411
00412 void NodalDOFMapHN::getDOFsForHNCell(
00413 int cellDim,
00414 int cellLID,
00415 int funcID,
00416 Array<int>& dofs ,
00417 Array<double>& coefs ) const {
00418
00419
00420 if ( (cellDim == 0) && (hangingNodeLID_to_NodesLIDs_.containsKey(cellLID)) )
00421 {
00422 Array<int> pointLIDs;
00423 Array<int> pointCoefs;
00424 pointLIDs = hangingNodeLID_to_NodesLIDs_.get( cellLID );
00425 dofs.resize(pointLIDs.size());
00426
00427 for (int ind = 0 ; ind < pointLIDs.size() ; ind++){
00428 dofs[ind] = nodeDofs_[ pointLIDs[ind] *nFuncs_ + funcID ];
00429 }
00430
00431 coefs = hangindNodeLID_to_Coefs_.get( cellLID );
00432 }
00433 }
00434
00435
00436
00437 void NodalDOFMapHN::getPointLIDsForHN( int pointLID , int facetIndex ,
00438 int maxCellIndex ,Array<int>& glbLIDs, Array<double>& coefsArray , Array<int>& nodeIndex){
00439
00440 Array<int> facets;
00441 int indexInParent , parentLID , facetCase = 0;
00442 double divRatio = 0.5;
00443 switch (dim_){
00444 case 2:
00445
00446 glbLIDs.resize(2);
00447 nodeIndex.resize(2);
00448 indexInParent = mesh().indexInParent(maxCellIndex);
00449 switch (indexInParent){
00450 case 0: {facetCase = (facetIndex == 1)? 0 : 1;
00451 divRatio = (1.0/3.0); break;}
00452 case 1: {facetCase = 0;
00453 divRatio = (facetIndex == 0)? (1.0/3.0) : (2.0/3.0); break;}
00454 case 2: {facetCase = (facetIndex == 0)? 0 : 2;
00455 divRatio = (facetIndex == 0)? (2.0/3.0) : (1.0/3.0); break;}
00456 case 5: {facetCase = 1;
00457 divRatio = (facetIndex == 0)? (1.0/3.0) : (2.0/3.0); break;}
00458 case 3: {facetCase = 2;
00459 divRatio = (facetIndex == 1)? (1.0/3.0) : (2.0/3.0); break;}
00460 case 6: {facetCase = (facetIndex == 0)? 1 : 3;
00461 divRatio = (facetIndex == 0)? (2.0/3.0) : (1.0/3.0); break;}
00462 case 7: {facetCase = 3;
00463 divRatio = (facetIndex == 2)? (1.0/3.0) : (2.0/3.0); break;}
00464 case 8: {facetCase = (facetIndex == 1)? 2 : 3;
00465 divRatio = (2.0/3.0); break;}
00466 }
00467 mesh().returnParentFacets(maxCellIndex , 0 ,facets , parentLID );
00468 switch (facetCase){
00469 case 0: {glbLIDs[0] = facets[0]; glbLIDs[1] = facets[1];
00470 nodeIndex[0] = 0; nodeIndex[1] = 1; break;}
00471 case 1: {glbLIDs[0] = facets[0]; glbLIDs[1] = facets[2];
00472 nodeIndex[0] = 0; nodeIndex[1] = 2; break;}
00473 case 2: {glbLIDs[0] = facets[1]; glbLIDs[1] = facets[3];
00474 nodeIndex[0] = 1; nodeIndex[1] = 3; break;}
00475 case 3: {glbLIDs[0] = facets[2]; glbLIDs[1] = facets[3];
00476 nodeIndex[0] = 2; nodeIndex[1] = 3; break;}
00477 }
00478 coefsArray.resize(2); coefsArray[0] = 1 - divRatio; coefsArray[1] = divRatio;
00479 SUNDANCE_MSG2(setupVerb(),"NodalDOFMapHN::getPointLIDsForHN() fc=" << facetCase << " R=" << divRatio
00480 << " facetIndex:" << facetIndex << " indexInParent:" << indexInParent <<" glbLIDs:"
00481 << glbLIDs << " maxCellIndex:" << maxCellIndex << " nodeIndex:" << nodeIndex);
00482 break;
00483 case 3:{
00484
00485 double ind_x = 0 , ind_y = 0 , ind_z = 0;
00486 double f_x = 0.0 , f_y = 0.0 , f_z = 0.0;
00487 double xx = 0.0 , yy = 0.0 , zz = 0.0;
00488 double values[8];
00489 indexInParent = mesh().indexInParent(maxCellIndex);
00490 mesh().returnParentFacets(maxCellIndex , 0 ,facets , parentLID );
00491
00492 switch (indexInParent){
00493 case 0: { ind_x = 0.0; ind_y = 0.0; ind_z = 0.0; break;}
00494 case 1: { ind_x = 1.0; ind_y = 0.0; ind_z = 0.0; break;}
00495 case 2: { ind_x = 2.0; ind_y = 0.0; ind_z = 0.0; break;}
00496 case 3: { ind_x = 0.0; ind_y = 1.0; ind_z = 0.0; break;}
00497 case 4: { ind_x = 1.0; ind_y = 1.0; ind_z = 0.0; break;}
00498 case 5: { ind_x = 2.0; ind_y = 1.0; ind_z = 0.0; break;}
00499 case 6: { ind_x = 0.0; ind_y = 2.0; ind_z = 0.0; break;}
00500 case 7: { ind_x = 1.0; ind_y = 2.0; ind_z = 0.0; break;}
00501 case 8: { ind_x = 2.0; ind_y = 2.0; ind_z = 0.0; break;}
00502 case 9: { ind_x = 0.0; ind_y = 0.0; ind_z = 1.0; break;}
00503 case 10: { ind_x = 1.0; ind_y = 0.0; ind_z = 1.0; break;}
00504 case 11: { ind_x = 2.0; ind_y = 0.0; ind_z = 1.0; break;}
00505 case 12: { ind_x = 0.0; ind_y = 1.0; ind_z = 1.0; break;}
00506 case 14: { ind_x = 2.0; ind_y = 1.0; ind_z = 1.0; break;}
00507 case 15: { ind_x = 0.0; ind_y = 2.0; ind_z = 1.0; break;}
00508 case 16: { ind_x = 1.0; ind_y = 2.0; ind_z = 1.0; break;}
00509 case 17: { ind_x = 2.0; ind_y = 2.0; ind_z = 1.0; break;}
00510 case 18: { ind_x = 0.0; ind_y = 0.0; ind_z = 2.0; break;}
00511 case 19: { ind_x = 1.0; ind_y = 0.0; ind_z = 2.0; break;}
00512 case 20: { ind_x = 2.0; ind_y = 0.0; ind_z = 2.0; break;}
00513 case 21: { ind_x = 0.0; ind_y = 1.0; ind_z = 2.0; break;}
00514 case 22: { ind_x = 1.0; ind_y = 1.0; ind_z = 2.0; break;}
00515 case 23: { ind_x = 2.0; ind_y = 1.0; ind_z = 2.0; break;}
00516 case 24: { ind_x = 0.0; ind_y = 2.0; ind_z = 2.0; break;}
00517 case 25: { ind_x = 1.0; ind_y = 2.0; ind_z = 2.0; break;}
00518 case 26: { ind_x = 2.0; ind_y = 2.0; ind_z = 2.0; break;}
00519 default: {}
00520 }
00521 switch (facetIndex){
00522 case 0: { f_x = 0.0 , f_y = 0.0 , f_z = 0.0; break;}
00523 case 1: { f_x = 1.0 , f_y = 0.0 , f_z = 0.0; break;}
00524 case 2: { f_x = 0.0 , f_y = 1.0 , f_z = 0.0; break;}
00525 case 3: { f_x = 1.0 , f_y = 1.0 , f_z = 0.0; break;}
00526 case 4: { f_x = 0.0 , f_y = 0.0 , f_z = 1.0; break;}
00527 case 5: { f_x = 1.0 , f_y = 0.0 , f_z = 1.0; break;}
00528 case 6: { f_x = 0.0 , f_y = 1.0 , f_z = 1.0; break;}
00529 case 7: { f_x = 1.0 , f_y = 1.0 , f_z = 1.0; break;}
00530 default: {}
00531 }
00532
00533 xx = (ind_x + f_x)/3.0;
00534 yy = (ind_y + f_y)/3.0;
00535 zz = (ind_z + f_z)/3.0;
00536 values[0] = (1.0 - xx)*(1.0 - yy)*(1.0 - zz);
00537 values[1] = (xx)*(1.0 - yy)*(1.0 - zz);
00538 values[2] = (1.0 - xx)*(yy)*(1.0 - zz);
00539 values[3] = (xx)*(yy)*(1.0 - zz);
00540 values[4] = (1.0 - xx)*(1.0 - yy)*(zz);
00541 values[5] = (xx)*(1.0 - yy)*(zz);
00542 values[6] = (1.0 - xx)*(yy)*(zz);
00543 values[7] = (xx)*(yy)*(zz);
00544
00545 glbLIDs.resize(0);
00546 nodeIndex.resize(0);
00547 coefsArray.resize(0);
00548 for (int ii = 0 ; ii < 8 ; ii++ ){
00549
00550 if (fabs(values[ii]) > 1e-5){
00551 glbLIDs.append(facets[ii]);
00552 nodeIndex.append(ii);
00553 coefsArray.append(values[ii]);
00554 }
00555 }
00556 SUNDANCE_MSG2(setupVerb(),"NodalDOFMapHN::getPointLIDsForHN()" << " facetIndex:" << facetIndex << " indexInParent:"
00557 << indexInParent <<" glbLIDs:" << glbLIDs << " maxCellIndex:" << maxCellIndex << " nodeIndex:" << nodeIndex);
00558 break;
00559 }
00560 }
00561
00562 }
00563
00564
00565
00566 void NodalDOFMapHN::computeOffsets(int localCount)
00567 {
00568 Array<int> dofOffsets;
00569 int totalDOFCount;
00570 int myOffset = 0;
00571 int np = mesh().comm().getNProc();
00572 if (np > 1)
00573 {
00574 SUNDANCE_MSG2(setupVerb(),"NodalDOFMapHN::computeOffsets, localCount:" << localCount);
00575 MPIContainerComm<int>::accumulate(localCount, dofOffsets, totalDOFCount,
00576 mesh().comm());
00577 myOffset = dofOffsets[mesh().comm().getRank()];
00578
00579 SUNDANCE_MSG2(setupVerb(),"NodalDOFMapHN::computeOffsets, offset:" << myOffset);
00580 int nDofs = nNodes_ * nFuncs_;
00581 for (int i=0; i<nDofs; i++)
00582 {
00583 if (nodeDofs_[i] >= 0) nodeDofs_[i] += myOffset;
00584 }
00585 }
00586 else
00587 {
00588 totalDOFCount = localCount;
00589 }
00590
00591 setLowestLocalDOF(myOffset);
00592 setNumLocalDOFs(localCount);
00593 setTotalNumDOFs(totalDOFCount);
00594 }
00595
00596
00597 void NodalDOFMapHN::shareRemoteDOFs(const Array<Array<int> >& outgoingCellRequests)
00598 {
00599 int np = mesh().comm().getNProc();
00600 if (np==1) return;
00601 int rank = mesh().comm().getRank();
00602
00603 Array<Array<int> > incomingCellRequests;
00604 Array<Array<int> > outgoingDOFs(np);
00605 Array<Array<int> > incomingDOFs;
00606
00607 SUNDANCE_MSG2(setupVerb(),
00608 "p=" << mesh().comm().getRank()
00609 << "synchronizing DOFs for cells of dimension 0");
00610 SUNDANCE_MSG2(setupVerb(),
00611 "p=" << mesh().comm().getRank()
00612 << " sending cell reqs d=0, GID="
00613 << outgoingCellRequests);
00614
00615
00616 MPIContainerComm<int>::allToAll(outgoingCellRequests,
00617 incomingCellRequests,
00618 mesh().comm());
00619
00620
00621
00622 for (int p=0; p<np; p++)
00623 {
00624 if (p==rank) continue;
00625 const Array<int>& requestsFromProc = incomingCellRequests[p];
00626 int nReq = requestsFromProc.size();
00627
00628 SUNDANCE_MSG3(setupVerb(),"p=" << mesh().comm().getRank()
00629 << " recv'd from proc=" << p
00630 << " reqs for DOFs for cells "
00631 << requestsFromProc);
00632
00633 outgoingDOFs[p].resize(nReq);
00634
00635 for (int c=0; c<nReq; c++)
00636 {
00637 int GID = requestsFromProc[c];
00638 SUNDANCE_MSG2(setupVerb(),
00639 "p=" << rank
00640 << " processing zero-cell with GID=" << GID);
00641 int LID = mesh().mapGIDToLID(0, GID);
00642 SUNDANCE_MSG2(setupVerb(),
00643 "p=" << rank
00644 << " LID=" << LID << " dofs="
00645 << nodeDofs_[LID*nFuncs_]);
00646 outgoingDOFs[p][c] = nodeDofs_[LID*nFuncs_];
00647 SUNDANCE_MSG2(setupVerb(),
00648 "p=" << rank
00649 << " done processing cell with GID=" << GID);
00650 }
00651 }
00652
00653 SUNDANCE_MSG2(setupVerb(),
00654 "p=" << mesh().comm().getRank()
00655 << "answering DOF requests for cells of dimension 0");
00656
00657
00658 MPIContainerComm<int>::allToAll(outgoingDOFs,
00659 incomingDOFs,
00660 mesh().comm());
00661
00662 SUNDANCE_MSG2(setupVerb(),
00663 "p=" << mesh().comm().getRank()
00664 << "communicated DOF answers for cells of dimension 0" );
00665
00666
00667
00668
00669 for (int p=0; p<mesh().comm().getNProc(); p++)
00670 {
00671 if (p==mesh().comm().getRank()) continue;
00672 const Array<int>& dofsFromProc = incomingDOFs[p];
00673 int numCells = dofsFromProc.size();
00674 for (int c=0; c<numCells; c++)
00675 {
00676 int cellGID = outgoingCellRequests[p][c];
00677 int cellLID = mesh().mapGIDToLID(0, cellGID);
00678 int dof = dofsFromProc[c];
00679 for (int i=0; i<nFuncs_; i++)
00680 {
00681 nodeDofs_[cellLID*nFuncs_ + i] = dof+i;
00682 addGhostIndex(dof+i);
00683 }
00684 }
00685 }
00686 }