MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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  const bool coupled,
64  const int NumDimensions,
65  const int interpolationOrder,
66  const Array<GO> GFineNodesPerDir,
67  const Array<LO> LFineNodesPerDir) :
68  comm_(comm), coupled_(coupled), numDimensions(NumDimensions),
69  interpolationOrder_(interpolationOrder), gFineNodesPerDir(GFineNodesPerDir),
70  lFineNodesPerDir(LFineNodesPerDir) {
71 
72  coarseRate.resize(3);
73  endRate.resize(3);
77 
78  offsets.resize(3);
82 
83  } // Constructor
84 
85  template<class LocalOrdinal, class GlobalOrdinal, class Node>
87 
89  if(const char* dbg = std::getenv("MUELU_INDEXMANAGER_DEBUG")) {
90  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
91  out->setShowAllFrontMatter(false).setShowProcRank(true);
92  } else {
93  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
94  }
95 
96  if(coupled_) {
97  gNumFineNodes10 = gFineNodesPerDir[1]*gFineNodesPerDir[0];
98  gNumFineNodes = gFineNodesPerDir[2]*gNumFineNodes10;
99  } else {
100  gNumFineNodes10 = Teuchos::OrdinalTraits<GO>::invalid();
101  gNumFineNodes = Teuchos::OrdinalTraits<GO>::invalid();
102  }
103  lNumFineNodes10 = lFineNodesPerDir[1]*lFineNodesPerDir[0];
104  lNumFineNodes = lFineNodesPerDir[2]*lNumFineNodes10;
105  for(int dim = 0; dim < 3; ++dim) {
106  if(dim < numDimensions) {
107  if(coupled_) {
108  if(startIndices[dim] == 0) {
109  meshEdge[2*dim] = true;
110  }
111  if(startIndices[dim + 3] + 1 == gFineNodesPerDir[dim]) {
112  meshEdge[2*dim + 1] = true;
113  endRate[dim] = startIndices[dim + 3] % coarseRate[dim];
114  }
115  } else { // With uncoupled problem each rank might require a different endRate
116  meshEdge[2*dim] = true;
117  meshEdge[2*dim + 1] = true;
118  endRate[dim] = (lFineNodesPerDir[dim] - 1) % coarseRate[dim];
119  }
120  if(endRate[dim] == 0) {endRate[dim] = coarseRate[dim];}
121 
122  // If uncoupled aggregation is used, offsets[dim] = 0, so nothing to do.
123  if(coupled_) {
124  offsets[dim] = Teuchos::as<LO>(startIndices[dim]) % coarseRate[dim];
125  if(offsets[dim] == 0) {
126  coarseNodeOffsets[dim] = 0;
127  } else if(startIndices[dim] + endRate[dim] == lFineNodesPerDir[dim]) {
128  coarseNodeOffsets[dim] = endRate[dim] - offsets[dim];
129  } else {
130  coarseNodeOffsets[dim] = coarseRate[dim] - offsets[dim];
131  }
132 
133  if(interpolationOrder_ == 0) {
134  int rem = startIndices[dim] % coarseRate[dim];
135  if( (rem != 0) && (rem <= Teuchos::as<double>(coarseRate[dim]) / 2.0)) {
136  ghostInterface[2*dim] = true;
137  }
138  rem = startIndices[dim + 3] % coarseRate[dim];
139  // uncoupled by nature does not require ghosts nodes
140  if(coupled_ && (startIndices[dim + 3] != gFineNodesPerDir[dim] - 1) &&
141  (rem > Teuchos::as<double>(coarseRate[dim]) / 2.0)) {
142  ghostInterface[2*dim + 1] = true;
143  }
144 
145  } else if(interpolationOrder_ == 1) {
146  if(coupled_ && (startIndices[dim] % coarseRate[dim] != 0 ||
147  startIndices[dim] == gFineNodesPerDir[dim]-1)) {
148  ghostInterface[2*dim] = true;
149  }
150  if(coupled_ && (startIndices[dim + 3] != gFineNodesPerDir[dim] - 1) &&
151  ((lFineNodesPerDir[dim] == 1) || (startIndices[dim + 3] % coarseRate[dim] != 0))) {
152  ghostInterface[2*dim+1] = true;
153  }
154  }
155  }
156  } else { // Default value for dim >= numDimensions
157  endRate[dim] = 1;
158  }
159  }
160 
161  *out << "gFineNodesPerDir: " << gFineNodesPerDir << std::endl;
162  *out << "lFineNodesPerDir: " << lFineNodesPerDir << std::endl;
163  *out << "endRate: " << endRate << std::endl;
164  *out << "ghostInterface: {" << ghostInterface[0] << ", " << ghostInterface[1] << ", "
165  << ghostInterface[2] << ", " << ghostInterface[3] << ", " << ghostInterface[4] << ", "
166  << ghostInterface[5] << "}" << std::endl;
167  *out << "meshEdge: {" << meshEdge[0] << ", " << meshEdge[1] << ", "
168  << meshEdge[2] << ", " << meshEdge[3] << ", " << meshEdge[4] << ", "
169  << meshEdge[5] << "}" << std::endl;
170  *out << "startIndices: " << startIndices << std::endl;
171  *out << "offsets: " << offsets << std::endl;
172  *out << "coarseNodeOffsets: " << coarseNodeOffsets << std::endl;
173 
174  // Here one element can represent either the degenerate case of one node or the more general
175  // case of two nodes, i.e. x---x is a 1D element with two nodes and x is a 1D element with
176  // one node. This helps generating a 3D space from tensorial products...
177  // A good way to handle this would be to generalize the algorithm to take into account the
178  // discretization order used in each direction, at least in the FEM sense, since a 0 degree
179  // discretization will have a unique node per element. This way 1D discretization can be
180  // viewed as a 3D problem with one 0 degree element in the y direction and one 0 degre
181  // element in the z direction.
182  // !!! Operations below are aftecting both local and global values that have two !!!
183  // different orientations. Orientations can be interchanged using mapDirG2L and mapDirL2G.
184  // coarseRate, endRate and offsets are in the global basis, as well as all the variables
185  // starting with a g.
186  // !!! while the variables starting with an l are in the local basis. !!!
187  for(int dim = 0; dim < 3; ++dim) {
188  if(dim < numDimensions) {
189  // Check whether the partition includes the "end" of the mesh which means that endRate
190  // will apply. Also make sure that endRate is not 0 which means that the mesh does not
191  // require a particular treatment at the boundaries.
192  if( meshEdge[2*dim + 1] ) {
193  lCoarseNodesPerDir[dim] = (lFineNodesPerDir[dim] - endRate[dim] + offsets[dim] - 1)
194  / coarseRate[dim] + 1;
195  if(offsets[dim] == 0) {++lCoarseNodesPerDir[dim];}
196  } else {
197  lCoarseNodesPerDir[dim] = (lFineNodesPerDir[dim] + offsets[dim] - 1) / coarseRate[dim];
198  if(offsets[dim] == 0) {++lCoarseNodesPerDir[dim];}
199  }
200 
201  // The first branch of this if-statement will be used if the rank contains only one layer
202  // of nodes in direction i, that layer must also coincide with the boundary of the mesh
203  // and coarseRate[i] == endRate[i]...
204  if(interpolationOrder_ == 0) {
205  startGhostedCoarseNode[dim] = startIndices[dim] / coarseRate[dim];
206  int rem = startIndices[dim] % coarseRate[dim];
207  if(rem > (Teuchos::as<double>(coarseRate[dim]) / 2.0) ) {
208  ++startGhostedCoarseNode[dim];
209  }
210  } else {
211  if((startIndices[dim] == gFineNodesPerDir[dim] - 1) &&
212  (startIndices[dim] % coarseRate[dim] == 0)) {
213  startGhostedCoarseNode[dim] = startIndices[dim] / coarseRate[dim] - 1;
214  } else {
215  startGhostedCoarseNode[dim] = startIndices[dim] / coarseRate[dim];
216  }
217  }
218 
219  // This array is passed to the RAPFactory and eventually becomes gFineNodePerDir on the next
220  // level.
221  gCoarseNodesPerDir[dim] = (gFineNodesPerDir[dim] - 1) / coarseRate[dim];
222  if((gFineNodesPerDir[dim] - 1) % coarseRate[dim] == 0) {
223  ++gCoarseNodesPerDir[dim];
224  } else {
225  gCoarseNodesPerDir[dim] += 2;
226  }
227  } else { // Default value for dim >= numDimensions
228  // endRate[dim] = 1;
229  gCoarseNodesPerDir[dim] = 1;
230  lCoarseNodesPerDir[dim] = 1;
231  } // if (dim < numDimensions)
232 
233  // This would happen if the rank does not own any nodes but in that case a subcommunicator
234  // should be used so this should really not be a concern.
235  if(lFineNodesPerDir[dim] < 1) {lCoarseNodesPerDir[dim] = 0;}
236  ghostedNodesPerDir[dim] = lCoarseNodesPerDir[dim];
237  // Check whether face *low needs ghost nodes
238  if(ghostInterface[2*dim]) {ghostedNodesPerDir[dim] += 1;}
239  // Check whether face *hi needs ghost nodes
240  if(ghostInterface[2*dim + 1]) {ghostedNodesPerDir[dim] += 1;}
241  } // Loop for dim=0:3
242 
243  // With uncoupled aggregation we need to communicate to compute the global number of coarse points
244  if(!coupled_) {
245  for(int dim = 0; dim < 3; ++dim) {
246  gCoarseNodesPerDir[dim] = -1;
247  }
248  }
249 
250  // Compute cummulative values
251  lNumCoarseNodes10 = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[1];
252  lNumCoarseNodes = lNumCoarseNodes10*lCoarseNodesPerDir[2];
253  numGhostedNodes10 = ghostedNodesPerDir[1]*ghostedNodesPerDir[0];
254  numGhostedNodes = numGhostedNodes10*ghostedNodesPerDir[2];
255  numGhostNodes = numGhostedNodes - lNumCoarseNodes;
256 
257  *out << "lCoarseNodesPerDir: " << lCoarseNodesPerDir << std::endl;
258  *out << "gCoarseNodesPerDir: " << gCoarseNodesPerDir << std::endl;
259  *out << "ghostedNodesPerDir: " << ghostedNodesPerDir << std::endl;
260  *out << "gNumCoarseNodes=" << gNumCoarseNodes << std::endl;
261  *out << "lNumCoarseNodes=" << lNumCoarseNodes << std::endl;
262  *out << "numGhostedNodes=" << numGhostedNodes << std::endl;
263  }
264 
265 } //namespace MueLu
266 
267 #define MUELU_INDEXMANAGER_SHORT
268 #endif // MUELU_INDEXMANAGER_DEF_HPP
Array< GO > gCoarseNodesPerDir
global number of nodes per direction remaining after coarsening.
basic_FancyOStream & setShowProcRank(const bool showProcRank)
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.
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
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
void resize(size_type new_size, const value_type &x=value_type())
Array< LO > coarseNodeOffsets
distance between lowest (resp. highest) index to the lowest (resp. highest) coarseNodeIndex in that d...