MueLu  Version of the Day
MueLu_IndexManager_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Ray Tuminaro (rstumin@sandia.gov)
41 // Luc Berger-Vergiat (lberge@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_INDEXMANAGER_DEF_HPP
47 #define MUELU_INDEXMANAGER_DEF_HPP
48 
49 #include "Teuchos_OrdinalTraits.hpp"
50 
51 #include "MueLu_ConfigDefs.hpp"
52 #include "MueLu_BaseClass.hpp"
54 
55 /*****************************************************************************
56 
57 ****************************************************************************/
58 
59 namespace MueLu {
60 
61  template<class LocalOrdinal, class GlobalOrdinal, class Node>
63  IndexManager(const RCP<const Teuchos::Comm<int> > comm,
64  const bool coupled,
65  const bool singleCoarsePoint,
66  const int NumDimensions,
67  const int interpolationOrder,
68  const Array<GO> GFineNodesPerDir,
69  const Array<LO> LFineNodesPerDir) :
70  comm_(comm), coupled_(coupled), singleCoarsePoint_(singleCoarsePoint),
71  numDimensions(NumDimensions), interpolationOrder_(interpolationOrder),
72  gFineNodesPerDir(GFineNodesPerDir), lFineNodesPerDir(LFineNodesPerDir) {
73 
74  coarseRate.resize(3);
75  endRate.resize(3);
76  gCoarseNodesPerDir.resize(3);
77  lCoarseNodesPerDir.resize(3);
78  ghostedNodesPerDir.resize(3);
79 
80  offsets.resize(3);
81  coarseNodeOffsets.resize(3);
82  startIndices.resize(6);
83  startGhostedCoarseNode.resize(3);
84 
85  } // Constructor
86 
87  template<class LocalOrdinal, class GlobalOrdinal, class Node>
90 
91  RCP<Teuchos::FancyOStream> out;
92  if(const char* dbg = std::getenv("MUELU_INDEXMANAGER_DEBUG")) {
93  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
94  out->setShowAllFrontMatter(false).setShowProcRank(true);
95  } else {
96  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
97  }
98 
99  if(coupled_) {
100  gNumFineNodes10 = gFineNodesPerDir[1]*gFineNodesPerDir[0];
101  gNumFineNodes = gFineNodesPerDir[2]*gNumFineNodes10;
102  } else {
103  gNumFineNodes10 = Teuchos::OrdinalTraits<GO>::invalid();
104  gNumFineNodes = Teuchos::OrdinalTraits<GO>::invalid();
105  }
106  lNumFineNodes10 = lFineNodesPerDir[1]*lFineNodesPerDir[0];
107  lNumFineNodes = lFineNodesPerDir[2]*lNumFineNodes10;
108  for(int dim = 0; dim < 3; ++dim) {
109  if(dim < numDimensions) {
110  if(coupled_) {
111  if(startIndices[dim] == 0) {
112  meshEdge[2*dim] = true;
113  }
114  if(startIndices[dim + 3] + 1 == gFineNodesPerDir[dim]) {
115  meshEdge[2*dim + 1] = true;
116  endRate[dim] = startIndices[dim + 3] % coarseRate[dim];
117  }
118  } else { // With uncoupled problem each rank might require a different endRate
119  meshEdge[2*dim] = true;
120  meshEdge[2*dim + 1] = true;
121  endRate[dim] = (lFineNodesPerDir[dim] - 1) % coarseRate[dim];
122  }
123  if(endRate[dim] == 0) {endRate[dim] = coarseRate[dim];}
124 
125  // If uncoupled aggregation is used, offsets[dim] = 0, so nothing to do.
126  if(coupled_) {
127  offsets[dim] = Teuchos::as<LO>(startIndices[dim]) % coarseRate[dim];
128  if(offsets[dim] == 0) {
129  coarseNodeOffsets[dim] = 0;
130  } else if(startIndices[dim] + endRate[dim] == lFineNodesPerDir[dim]) {
131  coarseNodeOffsets[dim] = endRate[dim] - offsets[dim];
132  } else {
133  coarseNodeOffsets[dim] = coarseRate[dim] - offsets[dim];
134  }
135 
136  if(interpolationOrder_ == 0) {
137  int rem = startIndices[dim] % coarseRate[dim];
138  if( (rem != 0) && (rem <= Teuchos::as<double>(coarseRate[dim]) / 2.0)) {
139  ghostInterface[2*dim] = true;
140  }
141  rem = startIndices[dim + 3] % coarseRate[dim];
142  // uncoupled by nature does not require ghosts nodes
143  if(coupled_ && (startIndices[dim + 3] != gFineNodesPerDir[dim] - 1) &&
144  (rem > Teuchos::as<double>(coarseRate[dim]) / 2.0)) {
145  ghostInterface[2*dim + 1] = true;
146  }
147 
148  } else if(interpolationOrder_ == 1) {
149  if(coupled_ && (startIndices[dim] % coarseRate[dim] != 0 ||
150  startIndices[dim] == gFineNodesPerDir[dim]-1)) {
151  ghostInterface[2*dim] = true;
152  }
153  if(coupled_ && (startIndices[dim + 3] != gFineNodesPerDir[dim] - 1) &&
154  ((lFineNodesPerDir[dim] == 1) || (startIndices[dim + 3] % coarseRate[dim] != 0))) {
155  ghostInterface[2*dim+1] = true;
156  }
157  }
158  }
159  } else { // Default value for dim >= numDimensions
160  endRate[dim] = 1;
161  }
162  }
163 
164  *out << "singleCoarsePoint? " << singleCoarsePoint_ << std::endl;
165  *out << "gFineNodesPerDir: " << gFineNodesPerDir << std::endl;
166  *out << "lFineNodesPerDir: " << lFineNodesPerDir << std::endl;
167  *out << "endRate: " << endRate << std::endl;
168  *out << "ghostInterface: {" << ghostInterface[0] << ", " << ghostInterface[1] << ", "
169  << ghostInterface[2] << ", " << ghostInterface[3] << ", " << ghostInterface[4] << ", "
170  << ghostInterface[5] << "}" << std::endl;
171  *out << "meshEdge: {" << meshEdge[0] << ", " << meshEdge[1] << ", "
172  << meshEdge[2] << ", " << meshEdge[3] << ", " << meshEdge[4] << ", "
173  << meshEdge[5] << "}" << std::endl;
174  *out << "startIndices: " << startIndices << std::endl;
175  *out << "offsets: " << offsets << std::endl;
176  *out << "coarseNodeOffsets: " << coarseNodeOffsets << std::endl;
177 
178  // Here one element can represent either the degenerate case of one node or the more general
179  // case of two nodes, i.e. x---x is a 1D element with two nodes and x is a 1D element with
180  // one node. This helps generating a 3D space from tensorial products...
181  // A good way to handle this would be to generalize the algorithm to take into account the
182  // discretization order used in each direction, at least in the FEM sense, since a 0 degree
183  // discretization will have a unique node per element. This way 1D discretization can be
184  // viewed as a 3D problem with one 0 degree element in the y direction and one 0 degre
185  // element in the z direction.
186  // !!! Operations below are aftecting both local and global values that have two !!!
187  // different orientations. Orientations can be interchanged using mapDirG2L and mapDirL2G.
188  // coarseRate, endRate and offsets are in the global basis, as well as all the variables
189  // starting with a g.
190  // !!! while the variables starting with an l are in the local basis. !!!
191  for(int dim = 0; dim < 3; ++dim) {
192  if(dim < numDimensions) {
193  // Check whether the partition includes the "end" of the mesh which means that endRate
194  // will apply. Also make sure that endRate is not 0 which means that the mesh does not
195  // require a particular treatment at the boundaries.
196  if( meshEdge[2*dim + 1] ) {
197  lCoarseNodesPerDir[dim] = (lFineNodesPerDir[dim] - endRate[dim] + offsets[dim] - 1)
198  / coarseRate[dim] + 1;
199  if(offsets[dim] == 0) {++lCoarseNodesPerDir[dim];}
200  // We might want to coarsening the direction
201  // into a single layer if there are not enough
202  // points left to form two aggregates
203  if(singleCoarsePoint_ && lFineNodesPerDir[dim] - 1 < coarseRate[dim]) {
204  lCoarseNodesPerDir[dim] =1;
205  }
206  } else {
207  lCoarseNodesPerDir[dim] = (lFineNodesPerDir[dim] + offsets[dim] - 1) / coarseRate[dim];
208  if(offsets[dim] == 0) {++lCoarseNodesPerDir[dim];}
209  }
210 
211  // The first branch of this if-statement will be used if the rank contains only one layer
212  // of nodes in direction i, that layer must also coincide with the boundary of the mesh
213  // and coarseRate[i] == endRate[i]...
214  if(interpolationOrder_ == 0) {
215  startGhostedCoarseNode[dim] = startIndices[dim] / coarseRate[dim];
216  int rem = startIndices[dim] % coarseRate[dim];
217  if(rem > (Teuchos::as<double>(coarseRate[dim]) / 2.0) ) {
218  ++startGhostedCoarseNode[dim];
219  }
220  } else {
221  if((startIndices[dim] == gFineNodesPerDir[dim] - 1) &&
222  (startIndices[dim] % coarseRate[dim] == 0)) {
223  startGhostedCoarseNode[dim] = startIndices[dim] / coarseRate[dim] - 1;
224  } else {
225  startGhostedCoarseNode[dim] = startIndices[dim] / coarseRate[dim];
226  }
227  }
228 
229  // This array is passed to the RAPFactory and eventually becomes gFineNodePerDir on the next
230  // level.
231  gCoarseNodesPerDir[dim] = (gFineNodesPerDir[dim] - 1) / coarseRate[dim];
232  if((gFineNodesPerDir[dim] - 1) % coarseRate[dim] == 0) {
233  ++gCoarseNodesPerDir[dim];
234  } else {
235  gCoarseNodesPerDir[dim] += 2;
236  }
237  } else { // Default value for dim >= numDimensions
238  // endRate[dim] = 1;
239  gCoarseNodesPerDir[dim] = 1;
240  lCoarseNodesPerDir[dim] = 1;
241  } // if (dim < numDimensions)
242 
243  // This would happen if the rank does not own any nodes but in that case a subcommunicator
244  // should be used so this should really not be a concern.
245  if(lFineNodesPerDir[dim] < 1) {lCoarseNodesPerDir[dim] = 0;}
246  ghostedNodesPerDir[dim] = lCoarseNodesPerDir[dim];
247  // Check whether face *low needs ghost nodes
248  if(ghostInterface[2*dim]) {ghostedNodesPerDir[dim] += 1;}
249  // Check whether face *hi needs ghost nodes
250  if(ghostInterface[2*dim + 1]) {ghostedNodesPerDir[dim] += 1;}
251  } // Loop for dim=0:3
252 
253  // With uncoupled aggregation we need to communicate to compute the global number of coarse points
254  if(!coupled_) {
255  for(int dim = 0; dim < 3; ++dim) {
256  gCoarseNodesPerDir[dim] = -1;
257  }
258  }
259 
260  // Compute cummulative values
261  lNumCoarseNodes10 = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[1];
262  lNumCoarseNodes = lNumCoarseNodes10*lCoarseNodesPerDir[2];
263  numGhostedNodes10 = ghostedNodesPerDir[1]*ghostedNodesPerDir[0];
264  numGhostedNodes = numGhostedNodes10*ghostedNodesPerDir[2];
265  numGhostNodes = numGhostedNodes - lNumCoarseNodes;
266 
267  *out << "lCoarseNodesPerDir: " << lCoarseNodesPerDir << std::endl;
268  *out << "gCoarseNodesPerDir: " << gCoarseNodesPerDir << std::endl;
269  *out << "ghostedNodesPerDir: " << ghostedNodesPerDir << std::endl;
270  *out << "gNumCoarseNodes=" << gNumCoarseNodes << std::endl;
271  *out << "lNumCoarseNodes=" << lNumCoarseNodes << std::endl;
272  *out << "numGhostedNodes=" << numGhostedNodes << std::endl;
273  }
274 
275 } //namespace MueLu
276 
277 #define MUELU_INDEXMANAGER_SHORT
278 #endif // MUELU_INDEXMANAGER_DEF_HPP
Array< GO > gCoarseNodesPerDir
global number of nodes per direction remaining after coarsening.
Array< GO > startIndices
lowest global tuple (i,j,k) of a node on the local process
Array< GO > startGhostedCoarseNode
lowest coarse global tuple (i,j,k) of a node remaing on the local process after coarsening.
Namespace for MueLu classes and methods.
Array< LO > offsets
distance between lowest (resp. highest) index to the lowest (resp. highest) ghostedNodeIndex in that ...
Array< LO > ghostedNodesPerDir
local number of ghosted nodes (i.e. ghost + coarse nodes) per direction
Array< int > coarseRate
coarsening rate in each direction
IndexManager()=default
Array< LO > lCoarseNodesPerDir
local number of nodes per direction remaing after coarsening.
Array< int > endRate
adapted coarsening rate at the edge of the mesh in each direction.
Array< LO > coarseNodeOffsets
distance between lowest (resp. highest) index to the lowest (resp. highest) coarseNodeIndex in that d...