MueLu Version of the Day
Loading...
Searching...
No Matches
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
59namespace 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 << "lNumCoarseNodes=" << lNumCoarseNodes << std::endl;
271 *out << "numGhostedNodes=" << numGhostedNodes << std::endl;
272 }
273
274} //namespace MueLu
275
276#define MUELU_INDEXMANAGER_SHORT
277#endif // MUELU_INDEXMANAGER_DEF_HPP
Array< LO > offsets
distance between lowest (resp. highest) index to the lowest (resp. highest) ghostedNodeIndex in that ...
Array< LO > coarseNodeOffsets
distance between lowest (resp. highest) index to the lowest (resp. highest) coarseNodeIndex in that d...
Array< LO > ghostedNodesPerDir
local number of ghosted nodes (i.e. ghost + coarse nodes) per direction
Array< GO > startGhostedCoarseNode
lowest coarse global tuple (i,j,k) of a node remaing on the local process after coarsening.
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< int > coarseRate
coarsening rate in each direction
IndexManager()=default
Array< int > endRate
adapted coarsening rate at the edge of the mesh in each direction.
Array< LO > lCoarseNodesPerDir
local number of nodes per direction remaing after coarsening.
Namespace for MueLu classes and methods.