MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_Q2Q1Q2CoarseGridFactory_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 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_Q2Q1Q2COARSEGRIDFACTORY_DEF_HPP
47 #define MUELU_Q2Q1Q2COARSEGRIDFACTORY_DEF_HPP
48 
49 #include <iostream>
50 #include <cmath>
51 
53 
54 #include <Xpetra_Map.hpp>
55 #include <Xpetra_MapFactory.hpp>
56 #include <Xpetra_Matrix.hpp>
57 #include <Xpetra_MultiVector.hpp>
58 #include <Xpetra_MultiVectorFactory.hpp>
59 #include <Xpetra_VectorFactory.hpp>
60 #include <Xpetra_CrsMatrixWrap.hpp>
61 
63 #include <MueLu_Level.hpp>
64 
65 namespace MueLu {
66 
67 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69  GetOStream(Runtime1) << "I constructed a Q2Q1Q2CoarseGridFactory object... Nothing else to do here." << std::endl;
70 }
71 
72 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74  // Should be empty. All destruction should be handled by Level-based get stuff and RCP
75 }
76 
77 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
79  Input(fineLevel, "VElementList");
80  Input(fineLevel, "PElementList");
81  Input(fineLevel, "MElementList");
82 
83  Input(coarseLevel, "VElementList");
84  Input(coarseLevel, "PElementList");
85  Input(coarseLevel, "MElementList");
86 
87  // currentLevel.DeclareInput(varName_,factory_,this);
88 }
89 
90 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
92  GetOStream(Runtime1) << "Starting 'build' routine...\n";
93 
94  // This will create a list of elements on the coarse grid with a
95  // predictable structure, as well as modify the fine grid list of
96  // elements, if necessary (i.e. if fineLevel.GetLevelID()==0);
97  // BuildCoarseGrid(fineLevel,coarseLevel);
98 
99  // This will actually build our prolongator P
100  return BuildCoarseGrid(fineLevel, coarseLevel);
101 }
102 
103 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
105  GetOStream(Runtime1) << "starting 'BuildCoarseGrid' routine...\n";
106 
107  RCP<Teuchos::SerialDenseMatrix<GO, GO> > fineElementPDOFs = Get<RCP<Teuchos::SerialDenseMatrix<GO, GO> > >(fineLevel, "PElementList");
108 
109  GO totalFineElements = fineElementPDOFs->numRows();
110 
111  // Compute number of coarse grid elements in total:
112  GO totalCoarseElements = totalFineElements / 4;
113  LO nCoarseElements = (int)sqrt(totalCoarseElements);
114 
115  // Initialize some counters:
116  size_t EdgeCount = (nCoarseElements + 1) * (nCoarseElements + 1);
117  size_t CenterCount = EdgeCount + 2 * nCoarseElements * (nCoarseElements + 1);
118 
119  // Initialize arrays of the proper size:
120  RCP<Teuchos::SerialDenseMatrix<GO, GO> > coarseElementVDOFs = rcp(new Teuchos::SerialDenseMatrix<GO, GO>(totalCoarseElements, 18));
121  RCP<Teuchos::SerialDenseMatrix<GO, GO> > coarseElementPDOFs = rcp(new Teuchos::SerialDenseMatrix<GO, GO>(totalCoarseElements, 4));
122  RCP<Teuchos::SerialDenseMatrix<GO, GO> > coarseElementMDOFs = rcp(new Teuchos::SerialDenseMatrix<GO, GO>(totalCoarseElements, 9));
123 
124  for (GO coarseElement = 0; coarseElement < totalCoarseElements; coarseElement++) {
125  // ***************************************************************
126  // This is less of a pain in the ass for magnetics, so I'm
127  // going to build the magnetics list. The velocity follows
128  // by doubling everything (and adding 1 for the y-components)
129  // and the pressure follows by copying the magnetics nodes.
130  // ***************************************************************
131 
132  // if (coarseElement is on the Bottom Edge)
133  if (coarseElement < nCoarseElements) {
134  // Bottom nodes
135  (*coarseElementMDOFs)(coarseElement, 0) = coarseElement;
136  (*coarseElementMDOFs)(coarseElement, 1) = coarseElement + 1;
137 
138  // Bottom edge
139  (*coarseElementMDOFs)(coarseElement, 4) = EdgeCount++;
140 
141  } else {
142  // Bottom Nodes
143  (*coarseElementMDOFs)(coarseElement, 0) = (*coarseElementMDOFs)(coarseElement - nCoarseElements, 3);
144  (*coarseElementMDOFs)(coarseElement, 1) = (*coarseElementMDOFs)(coarseElement - nCoarseElements, 2);
145 
146  // Bottom Edge
147  (*coarseElementMDOFs)(coarseElement, 4) = (*coarseElementMDOFs)(coarseElement - nCoarseElements, 6);
148  }
149 
150  // Right and Top Edges -- must be determined before left edge
151  (*coarseElementMDOFs)(coarseElement, 5) = EdgeCount++;
152  (*coarseElementMDOFs)(coarseElement, 6) = EdgeCount++;
153 
154  // if (coarseElement is on the Left Edge)
155  if (coarseElement % nCoarseElements == 0) {
156  // Top left node
157  (*coarseElementMDOFs)(coarseElement, 3) = (*coarseElementMDOFs)(coarseElement, 0) + nCoarseElements + 1;
158 
159  // Left Edge
160  (*coarseElementMDOFs)(coarseElement, 7) = EdgeCount++;
161 
162  } else {
163  // Top left node
164  (*coarseElementMDOFs)(coarseElement, 3) = (*coarseElementMDOFs)(coarseElement - 1, 2);
165 
166  // Left Edge
167  (*coarseElementMDOFs)(coarseElement, 7) = (*coarseElementMDOFs)(coarseElement - 1, 5);
168  }
169 
170  // Top right node -- Must be the last node to be determined!
171  (*coarseElementMDOFs)(coarseElement, 2) = (*coarseElementMDOFs)(coarseElement, 3) + 1;
172 
173  // Center Node
174  (*coarseElementMDOFs)(coarseElement, 8) = CenterCount++;
175 
176  // With Magnetics built, Pressure and Velocity follow without effort.
177  // First, Velocity:
178  (*coarseElementVDOFs)(coarseElement, 0) = 2 * (*coarseElementMDOFs)(coarseElement, 0);
179  (*coarseElementVDOFs)(coarseElement, 1) = 2 * (*coarseElementMDOFs)(coarseElement, 0) + 1;
180  (*coarseElementVDOFs)(coarseElement, 2) = 2 * (*coarseElementMDOFs)(coarseElement, 1);
181  (*coarseElementVDOFs)(coarseElement, 3) = 2 * (*coarseElementMDOFs)(coarseElement, 1) + 1;
182  (*coarseElementVDOFs)(coarseElement, 4) = 2 * (*coarseElementMDOFs)(coarseElement, 2);
183  (*coarseElementVDOFs)(coarseElement, 5) = 2 * (*coarseElementMDOFs)(coarseElement, 2) + 1;
184  (*coarseElementVDOFs)(coarseElement, 6) = 2 * (*coarseElementMDOFs)(coarseElement, 3);
185  (*coarseElementVDOFs)(coarseElement, 7) = 2 * (*coarseElementMDOFs)(coarseElement, 3) + 1;
186  (*coarseElementVDOFs)(coarseElement, 8) = 2 * (*coarseElementMDOFs)(coarseElement, 4);
187  (*coarseElementVDOFs)(coarseElement, 9) = 2 * (*coarseElementMDOFs)(coarseElement, 4) + 1;
188  (*coarseElementVDOFs)(coarseElement, 10) = 2 * (*coarseElementMDOFs)(coarseElement, 5);
189  (*coarseElementVDOFs)(coarseElement, 11) = 2 * (*coarseElementMDOFs)(coarseElement, 5) + 1;
190  (*coarseElementVDOFs)(coarseElement, 12) = 2 * (*coarseElementMDOFs)(coarseElement, 6);
191  (*coarseElementVDOFs)(coarseElement, 13) = 2 * (*coarseElementMDOFs)(coarseElement, 6) + 1;
192  (*coarseElementVDOFs)(coarseElement, 14) = 2 * (*coarseElementMDOFs)(coarseElement, 7);
193  (*coarseElementVDOFs)(coarseElement, 15) = 2 * (*coarseElementMDOFs)(coarseElement, 7) + 1;
194  (*coarseElementVDOFs)(coarseElement, 16) = 2 * (*coarseElementMDOFs)(coarseElement, 8);
195  (*coarseElementVDOFs)(coarseElement, 17) = 2 * (*coarseElementMDOFs)(coarseElement, 8) + 1;
196 
197  // Lastly, Pressure:
198  (*coarseElementPDOFs)(coarseElement, 0) = (*coarseElementMDOFs)(coarseElement, 0);
199  (*coarseElementPDOFs)(coarseElement, 1) = (*coarseElementMDOFs)(coarseElement, 1);
200  (*coarseElementPDOFs)(coarseElement, 2) = (*coarseElementMDOFs)(coarseElement, 2);
201  (*coarseElementPDOFs)(coarseElement, 3) = (*coarseElementMDOFs)(coarseElement, 3);
202 
203  } // Loop over elements
204 
205  Set(coarseLevel, "VElementList", coarseElementVDOFs);
206  Set(coarseLevel, "PElementList", coarseElementPDOFs);
207  Set(coarseLevel, "MElementList", coarseElementMDOFs);
208 
209  // coarseLevel.Keep("VElementList",coarseLevel.GetFactoryManager()->GetFactory("VElementList").get());
210  // coarseLevel.Keep("PElementList",coarseLevel.GetFactoryManager()->GetFactory("PElementList").get());
211  // coarseLevel.Keep("MElementList",coarseLevel.GetFactoryManager()->GetFactory("MElementList").get());
212 
213 } // BuildCoarseGrid
214 
215 } // namespace MueLu
216 
217 #define MUELU_Q2Q1Q2COARSEGRIDFACTORY_SHORT
218 #endif // MUELU_Q2Q1Q2COARSEGRIDFACTORY_DEF_HPP
GlobalOrdinal GO
LocalOrdinal LO
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Description of what is happening (more verbose)
void BuildCoarseGrid(Level &fineLevel, Level &coarseLevel) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.