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>
59 #include <Xpetra_VectorFactory.hpp>
60 #include <Xpetra_CrsMatrixWrap.hpp>
61 #include <Xpetra_MultiVector.hpp>
63 
65 #include <MueLu_Level.hpp>
66 #include <MueLu_Monitor.hpp>
67 #include <MueLu_Utilities.hpp>
70 
71 namespace MueLu {
72 
73  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
75  GetOStream(Runtime1) << "I constructed a Q2Q1Q2CoarseGridFactory object... Nothing else to do here." << std::endl;
76  }
77 
78  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
80  // Should be empty. All destruction should be handled by Level-based get stuff and RCP
81  }
82 
83  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
85 
86  Input(fineLevel, "VElementList");
87  Input(fineLevel, "PElementList");
88  Input(fineLevel, "MElementList");
89 
90  Input(coarseLevel, "VElementList");
91  Input(coarseLevel, "PElementList");
92  Input(coarseLevel, "MElementList");
93 
94 
95  //currentLevel.DeclareInput(varName_,factory_,this);
96  }
97 
98  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
100 
101  GetOStream(Runtime1) << "Starting 'build' routine...\n";
102 
103  // This will create a list of elements on the coarse grid with a
104  // predictable structure, as well as modify the fine grid list of
105  // elements, if necessary (i.e. if fineLevel.GetLevelID()==0);
106  //BuildCoarseGrid(fineLevel,coarseLevel);
107 
108  // This will actually build our prolongator P
109  return BuildCoarseGrid(fineLevel,coarseLevel);
110 
111  }
112 
113 
114  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
116  {
117 
118  GetOStream(Runtime1) << "starting 'BuildCoarseGrid' routine...\n";
119 
120  RCP<Teuchos::SerialDenseMatrix<GO,GO> > fineElementPDOFs = Get< RCP<Teuchos::SerialDenseMatrix<GO,GO> > >(fineLevel,"PElementList");
121 
122  GO totalFineElements = fineElementPDOFs->numRows();
123 
124  // Compute number of coarse grid elements in total:
125  GO totalCoarseElements = totalFineElements/4;
126  LO nCoarseElements = (int) sqrt(totalCoarseElements);
127 
128  // Initialize some counters:
129  size_t EdgeCount = (nCoarseElements + 1) * (nCoarseElements + 1);
130  size_t CenterCount = EdgeCount + 2 * nCoarseElements * (nCoarseElements + 1);
131 
132 
133  // Initialize arrays of the proper size:
134  RCP<Teuchos::SerialDenseMatrix<GO,GO> > coarseElementVDOFs = rcp(new Teuchos::SerialDenseMatrix<GO,GO>(totalCoarseElements,18));
135  RCP<Teuchos::SerialDenseMatrix<GO,GO> > coarseElementPDOFs = rcp(new Teuchos::SerialDenseMatrix<GO,GO>(totalCoarseElements,4));
136  RCP<Teuchos::SerialDenseMatrix<GO,GO> > coarseElementMDOFs = rcp(new Teuchos::SerialDenseMatrix<GO,GO>(totalCoarseElements,9));
137 
138 
139  for ( GO coarseElement=0; coarseElement < totalCoarseElements; coarseElement++ )
140  {
141 
142  // ***************************************************************
143  // This is less of a pain in the ass for magnetics, so I'm
144  // going to build the magnetics list. The velocity follows
145  // by doubling everything (and adding 1 for the y-components)
146  // and the pressure follows by copying the magnetics nodes.
147  // ***************************************************************
148 
149  //if (coarseElement is on the Bottom Edge)
150  if (coarseElement < nCoarseElements)
151  {
152  //Bottom nodes
153  (*coarseElementMDOFs)(coarseElement,0) = coarseElement;
154  (*coarseElementMDOFs)(coarseElement,1) = coarseElement+1;
155 
156  //Bottom edge
157  (*coarseElementMDOFs)(coarseElement,4) = EdgeCount++;
158 
159  }
160  else
161  {
162  //Bottom Nodes
163  (*coarseElementMDOFs)(coarseElement,0) = (*coarseElementMDOFs)(coarseElement-nCoarseElements,3);
164  (*coarseElementMDOFs)(coarseElement,1) = (*coarseElementMDOFs)(coarseElement-nCoarseElements,2);
165 
166  //Bottom Edge
167  (*coarseElementMDOFs)(coarseElement,4) = (*coarseElementMDOFs)(coarseElement-nCoarseElements,6);
168 
169  }
170 
171 
172  //Right and Top Edges -- must be determined before left edge
173  (*coarseElementMDOFs)(coarseElement,5) = EdgeCount++;
174  (*coarseElementMDOFs)(coarseElement,6) = EdgeCount++;
175 
176 
177  //if (coarseElement is on the Left Edge)
178  if (coarseElement % nCoarseElements == 0)
179  {
180  //Top left node
181  (*coarseElementMDOFs)(coarseElement,3) = (*coarseElementMDOFs)(coarseElement,0)+nCoarseElements+1;
182 
183  //Left Edge
184  (*coarseElementMDOFs)(coarseElement,7) = EdgeCount++;
185 
186  }
187  else
188  {
189  //Top left node
190  (*coarseElementMDOFs)(coarseElement,3) = (*coarseElementMDOFs)(coarseElement-1,2);
191 
192  //Left Edge
193  (*coarseElementMDOFs)(coarseElement,7) = (*coarseElementMDOFs)(coarseElement-1,5);
194  }
195 
196  //Top right node -- Must be the last node to be determined!
197  (*coarseElementMDOFs)(coarseElement,2) = (*coarseElementMDOFs)(coarseElement,3)+1;
198 
199  //Center Node
200  (*coarseElementMDOFs)(coarseElement,8) = CenterCount++;
201 
202 
203 
204  // With Magnetics built, Pressure and Velocity follow without effort.
205  // First, Velocity:
206  (*coarseElementVDOFs)(coarseElement,0) = 2*(*coarseElementMDOFs)(coarseElement,0);
207  (*coarseElementVDOFs)(coarseElement,1) = 2*(*coarseElementMDOFs)(coarseElement,0)+1;
208  (*coarseElementVDOFs)(coarseElement,2) = 2*(*coarseElementMDOFs)(coarseElement,1);
209  (*coarseElementVDOFs)(coarseElement,3) = 2*(*coarseElementMDOFs)(coarseElement,1)+1;
210  (*coarseElementVDOFs)(coarseElement,4) = 2*(*coarseElementMDOFs)(coarseElement,2);
211  (*coarseElementVDOFs)(coarseElement,5) = 2*(*coarseElementMDOFs)(coarseElement,2)+1;
212  (*coarseElementVDOFs)(coarseElement,6) = 2*(*coarseElementMDOFs)(coarseElement,3);
213  (*coarseElementVDOFs)(coarseElement,7) = 2*(*coarseElementMDOFs)(coarseElement,3)+1;
214  (*coarseElementVDOFs)(coarseElement,8) = 2*(*coarseElementMDOFs)(coarseElement,4);
215  (*coarseElementVDOFs)(coarseElement,9) = 2*(*coarseElementMDOFs)(coarseElement,4)+1;
216  (*coarseElementVDOFs)(coarseElement,10) = 2*(*coarseElementMDOFs)(coarseElement,5);
217  (*coarseElementVDOFs)(coarseElement,11) = 2*(*coarseElementMDOFs)(coarseElement,5)+1;
218  (*coarseElementVDOFs)(coarseElement,12) = 2*(*coarseElementMDOFs)(coarseElement,6);
219  (*coarseElementVDOFs)(coarseElement,13) = 2*(*coarseElementMDOFs)(coarseElement,6)+1;
220  (*coarseElementVDOFs)(coarseElement,14) = 2*(*coarseElementMDOFs)(coarseElement,7);
221  (*coarseElementVDOFs)(coarseElement,15) = 2*(*coarseElementMDOFs)(coarseElement,7)+1;
222  (*coarseElementVDOFs)(coarseElement,16) = 2*(*coarseElementMDOFs)(coarseElement,8);
223  (*coarseElementVDOFs)(coarseElement,17) = 2*(*coarseElementMDOFs)(coarseElement,8)+1;
224 
225  // Lastly, Pressure:
226  (*coarseElementPDOFs)(coarseElement,0) = (*coarseElementMDOFs)(coarseElement,0);
227  (*coarseElementPDOFs)(coarseElement,1) = (*coarseElementMDOFs)(coarseElement,1);
228  (*coarseElementPDOFs)(coarseElement,2) = (*coarseElementMDOFs)(coarseElement,2);
229  (*coarseElementPDOFs)(coarseElement,3) = (*coarseElementMDOFs)(coarseElement,3);
230 
231  }// Loop over elements
232 
233  Set(coarseLevel,"VElementList",coarseElementVDOFs);
234  Set(coarseLevel,"PElementList",coarseElementPDOFs);
235  Set(coarseLevel,"MElementList",coarseElementMDOFs);
236 
237  //coarseLevel.Keep("VElementList",coarseLevel.GetFactoryManager()->GetFactory("VElementList").get());
238  //coarseLevel.Keep("PElementList",coarseLevel.GetFactoryManager()->GetFactory("PElementList").get());
239  //coarseLevel.Keep("MElementList",coarseLevel.GetFactoryManager()->GetFactory("MElementList").get());
240 
241  }//BuildCoarseGrid
242 
243 
244 } // namespace MueLu
245 
246 #define MUELU_Q2Q1Q2COARSEGRIDFACTORY_SHORT
247 #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.