MueLu Version of the Day
Loading...
Searching...
No Matches
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
52#include <Teuchos_SerialDenseMatrix.hpp>
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#include <Xpetra_MultiVector.hpp>
62#include <Xpetra_MultiVectorFactory.hpp>
63
65#include <MueLu_Level.hpp>
66#include <MueLu_Monitor.hpp>
67#include <MueLu_Utilities.hpp>
70
71namespace 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
Class that holds all level-specific information.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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.
Namespace for MueLu classes and methods.
@ Runtime1
Description of what is happening (more verbose)