MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_IntrepidPCoarsenFactory_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_IPCFACTORY_DEF_HPP
47#define MUELU_IPCFACTORY_DEF_HPP
48
49#include <Xpetra_Matrix.hpp>
50#include <Xpetra_IO.hpp>
51#include <sstream>
52#include <algorithm>
53
55
57#include "MueLu_Level.hpp"
58#include "MueLu_MasterList.hpp"
59#include "MueLu_Monitor.hpp"
60#include "MueLu_PerfUtils.hpp"
62#include "MueLu_Utilities.hpp"
63
64#include "Teuchos_ScalarTraits.hpp"
65
66// Intrepid Headers
67
68//Intrepid_HGRAD_HEX_C1_FEM.hpp
69//Intrepid_HGRAD_HEX_C2_FEM.hpp
70//Intrepid_HGRAD_HEX_Cn_FEM.hpp
71//Intrepid_HGRAD_HEX_I2_FEM.hpp
72#include "Intrepid2_HGRAD_LINE_C1_FEM.hpp"
73#include "Intrepid2_HGRAD_LINE_Cn_FEM.hpp"
74//Intrepid_HGRAD_LINE_Cn_FEM_JACOBI.hpp
75//Intrepid_HGRAD_POLY_C1_FEM.hpp
76//Intrepid_HGRAD_PYR_C1_FEM.hpp
77//Intrepid_HGRAD_PYR_I2_FEM.hpp
78#include "Intrepid2_HGRAD_QUAD_C1_FEM.hpp"
79//#include Intrepid_HGRAD_QUAD_C2_FEM.hpp
80#include "Intrepid2_HGRAD_QUAD_Cn_FEM.hpp"
81//Intrepid_HGRAD_TET_C1_FEM.hpp
82//Intrepid_HGRAD_TET_C2_FEM.hpp
83//Intrepid_HGRAD_TET_Cn_FEM.hpp
84//Intrepid_HGRAD_TET_Cn_FEM_ORTH.hpp
85//Intrepid_HGRAD_TET_COMP12_FEM.hpp
86//Intrepid_HGRAD_TRI_C1_FEM.hpp
87//Intrepid_HGRAD_TRI_C2_FEM.hpp
88//Intrepid_HGRAD_TRI_Cn_FEM.hpp
89//Intrepid_HGRAD_TRI_Cn_FEM_ORTH.hpp
90//Intrepid_HGRAD_WEDGE_C1_FEM.hpp
91//Intrepid_HGRAD_WEDGE_C2_FEM.hpp
92//Intrepid_HGRAD_WEDGE_I2_FEM.hpp
93
94// Helper Macro to avoid "unrequested" warnings
95#define MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(level,ename,entry) \
96 {if (level.IsRequested(ename,this) || level.GetKeepFlag(ename,this) != 0) this->Set(level,ename,entry);}
97
98
99
100namespace MueLu {
101
102
103/*********************************************************************************************************/
104namespace MueLuIntrepid {
105inline std::string tolower(const std::string & str) {
106 std::string data(str);
107 std::transform(data.begin(), data.end(), data.begin(), ::tolower);
108 return data;
109}
110
111
112/*********************************************************************************************************/
113 template<class Basis, class LOFieldContainer, class LocalOrdinal, class GlobalOrdinal, class Node>
114 void FindGeometricSeedOrdinals(Teuchos::RCP<Basis> basis, const LOFieldContainer &elementToNodeMap,
115 std::vector<std::vector<LocalOrdinal> > &seeds,
116 const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &rowMap,
117 const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &columnMap)
118 {
119
120 // For each subcell represented by the elements in elementToNodeMap, we want to identify a globally
121 // unique degree of freedom. Because the other "seed" interfaces in MueLu expect a local ordinal, we
122 // store local ordinals in the resulting seeds container.
123
124 // The approach is as follows. For each element, we iterate through the subcells of the domain topology.
125 // We determine which, if any, of these has the lowest global ID owned that is locally owned. We then insert
126 // the local ID corresponding to this in a vector<set<int>> container whose outer index is the spatial dimension
127 // of the subcell. The set lets us conveniently enforce uniqueness of the stored LIDs.
128
129 shards::CellTopology cellTopo = basis->getBaseCellTopology();
130 int spaceDim = cellTopo.getDimension();
131 seeds.clear();
132 seeds.resize(spaceDim + 1);
133 typedef GlobalOrdinal GO;
134 typedef LocalOrdinal LO;
135
136 LocalOrdinal lo_invalid = Teuchos::OrdinalTraits<LO>::invalid();
137 GlobalOrdinal go_invalid = Teuchos::OrdinalTraits<GO>::invalid();
138
139
140 std::vector<std::set<LocalOrdinal>> seedSets(spaceDim+1);
141
142 int numCells = elementToNodeMap.extent(0);
143 auto elementToNodeMap_host = Kokkos::create_mirror_view(elementToNodeMap);
144 Kokkos::deep_copy(elementToNodeMap_host,elementToNodeMap);
145 for (int cellOrdinal=0; cellOrdinal<numCells; cellOrdinal++)
146 {
147 for (int d=0; d<=spaceDim; d++)
148 {
149 int subcellCount = cellTopo.getSubcellCount(d);
150 for (int subcord=0; subcord<subcellCount; subcord++)
151 {
152 int dofCount = basis->getDofCount(d,subcord);
153 if (dofCount == 0) continue;
154 // otherwise, we want to insert the LID corresponding to the least globalID that is locally owned
155 GO leastGlobalDofOrdinal = go_invalid;
156 LO LID_leastGlobalDofOrdinal = lo_invalid;
157 for (int basisOrdinalOrdinal=0; basisOrdinalOrdinal<dofCount; basisOrdinalOrdinal++)
158 {
159 int basisOrdinal = basis->getDofOrdinal(d,subcord,basisOrdinalOrdinal);
160 int colLID = elementToNodeMap_host(cellOrdinal,basisOrdinal);
161 if (colLID != Teuchos::OrdinalTraits<LO>::invalid())
162 {
163 GlobalOrdinal colGID = columnMap.getGlobalElement(colLID);
164 LocalOrdinal rowLID = rowMap.getLocalElement(colGID);
165 if (rowLID != lo_invalid)
166 {
167 if ((leastGlobalDofOrdinal == go_invalid) || (colGID < leastGlobalDofOrdinal))
168 {
169 // replace with rowLID
170 leastGlobalDofOrdinal = colGID;
171 LID_leastGlobalDofOrdinal = rowLID;
172 }
173 }
174 }
175 }
176 if (leastGlobalDofOrdinal != go_invalid)
177 {
178 seedSets[d].insert(LID_leastGlobalDofOrdinal);
179 }
180 }
181 }
182 }
183 for (int d=0; d<=spaceDim;d++)
184 {
185 seeds[d] = std::vector<LocalOrdinal>(seedSets[d].begin(),seedSets[d].end());
186 }
187 }
188
189/*********************************************************************************************************/
190// Syntax [HGRAD|HCURL|HDIV][_| ][HEX|LINE|POLY|PYR|QUAD|TET|TRI|WEDGE][_| ][C|I][1|2|n]
191// Inputs:
192// name - name of the intrepid basis to generate
193// Outputs:
194// degree - order of resulting discretization
195// return value - Intrepid2 basis correspionding to the name
196template<class Scalar,class KokkosExecutionSpace>
197Teuchos::RCP< Intrepid2::Basis<KokkosExecutionSpace,Scalar,Scalar> > BasisFactory(const std::string & name, int & degree)
198{
199 using std::string;
200 using Teuchos::rcp;
201 string myerror("IntrepidBasisFactory: cannot parse string name '"+name+"'");
202
203 // Syntax [HGRAD|HCURL|HDIV][_| ][HEX|LINE|POLY|PYR|QUAD|TET|TRI|WEDGE][_| ][C|I][1|2|n]
204
205 // Get the derivative type
206 size_t pos1 = name.find_first_of(" _");
207 if(pos1==0) throw std::runtime_error(myerror);
208 string deriv = tolower(name.substr(0,pos1));
209 if(deriv!="hgrad" && deriv!="hcurl" && deriv!="hdiv") throw std::runtime_error(myerror);
210
211 // Get the element type
212 pos1++;
213 size_t pos2 = name.find_first_of(" _",pos1);
214 if(pos2==0) throw std::runtime_error(myerror);
215 string el = tolower(name.substr(pos1,pos2-pos1));
216 if(el!="hex" && el!="line" && el!="poly" && el!="pyr" && el!="quad" && el!="tet" && el!="tri" && el!="wedge") throw std::runtime_error(myerror);
217
218 // Get the polynomial type
219 pos2++;
220 string poly = tolower(name.substr(pos2,1));
221 if(poly!="c" && poly!="i") throw std::runtime_error(myerror);
222
223 // Get the degree
224 pos2++;
225 degree=std::stoi(name.substr(pos2,1));
226 if(degree<=0) throw std::runtime_error(myerror);
227
228 // FIXME LATER: Allow for alternative point types for Kirby elements
229 if(deriv=="hgrad" && el=="quad" && poly=="c"){
230 if(degree==1) return rcp(new Intrepid2::Basis_HGRAD_QUAD_C1_FEM<KokkosExecutionSpace,Scalar,Scalar>());
231 else return rcp(new Intrepid2::Basis_HGRAD_QUAD_Cn_FEM<KokkosExecutionSpace,Scalar,Scalar>(degree,Intrepid2::POINTTYPE_EQUISPACED));
232 }
233 else if(deriv=="hgrad" && el=="line" && poly=="c"){
234 if(degree==1) return rcp(new Intrepid2::Basis_HGRAD_LINE_C1_FEM<KokkosExecutionSpace,Scalar,Scalar>());
235 else return rcp(new Intrepid2::Basis_HGRAD_LINE_Cn_FEM<KokkosExecutionSpace,Scalar,Scalar>(degree,Intrepid2::POINTTYPE_EQUISPACED));
236 }
237
238 // Error out
239 throw std::runtime_error(myerror);
240 TEUCHOS_UNREACHABLE_RETURN(Teuchos::null);
241}
242
243/*********************************************************************************************************/
244// Gets the "lo" nodes nested into a "hi" basis. Only works on quads and lines for a lo basis of p=1
245// Inputs:
246// hi_basis - Higher order Basis
247// Outputs:
248// lo_node_in_hi - std::vector<size_t> of size lo dofs in the reference element, which describes the coindcident hi dots
249// hi_DofCoords - FC<Scalar> of size (#hi dofs, dim) with the coordinate locations of the hi dofs on the reference element
250template<class Scalar,class KokkosDeviceType>
251void IntrepidGetP1NodeInHi(const Teuchos::RCP<Intrepid2::Basis<typename KokkosDeviceType::execution_space,Scalar,Scalar> >&hi_basis,
252 std::vector<size_t> & lo_node_in_hi,
253 Kokkos::DynRankView<Scalar,KokkosDeviceType> & hi_DofCoords) {
254
255 typedef typename KokkosDeviceType::execution_space KokkosExecutionSpace;
256 // Figure out which unknowns in hi_basis correspond to nodes on lo_basis. This varies by element type.
257 size_t degree = hi_basis->getDegree();
258 lo_node_in_hi.resize(0);
259
260 if(!rcp_dynamic_cast<Intrepid2::Basis_HGRAD_QUAD_Cn_FEM<KokkosExecutionSpace,Scalar,Scalar> >(hi_basis).is_null()) {
261 // HGRAD QUAD Cn: Numbering as per the Kirby convention (straight across, bottom to top)
262 lo_node_in_hi.insert(lo_node_in_hi.end(),{0,degree, (degree+1)*(degree+1)-1, degree*(degree+1)});
263 }
264 else if(!rcp_dynamic_cast<Intrepid2::Basis_HGRAD_LINE_Cn_FEM<KokkosExecutionSpace,Scalar,Scalar> >(hi_basis).is_null()) {
265 // HGRAD LINE Cn: Numbering as per the Kirby convention (straight across)
266 lo_node_in_hi.insert(lo_node_in_hi.end(),{0,degree});
267 }
268 else
269 throw std::runtime_error("IntrepidPCoarsenFactory: Unknown element type");
270
271 // Get coordinates of the hi_basis dof's
272 Kokkos::resize(hi_DofCoords,hi_basis->getCardinality(),hi_basis->getBaseCellTopology().getDimension());
273 hi_basis->getDofCoords(hi_DofCoords);
274}
275
276
277/*********************************************************************************************************/
278// Given a list of candidates picks a definitive list of "representative" higher order nodes for each lo order node via the "smallest GID" rule
279// Input:
280// representative_node_candidates - std::vector<std::vector<size_t> > of lists of "representative candidate" hi dofs for each lo dof
281// hi_elemToNode - FC<LO> containing the high order element-to-node map
282// hi_columnMap - Column map of the higher order matrix
283// Output:
284// lo_elemToHiRepresentativeNode - FC<LO> of size (# elements, # lo dofs per element) listing the hi unknown chosen as the single representative for each lo unknown for counting purposes
285template<class LocalOrdinal, class GlobalOrdinal, class Node, class LOFieldContainer>
286void GenerateLoNodeInHiViaGIDs(const std::vector<std::vector<size_t> > & candidates,const LOFieldContainer & hi_elemToNode,
287 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > & hi_columnMap,
288 LOFieldContainer & lo_elemToHiRepresentativeNode) {
289 typedef GlobalOrdinal GO;
290
291 // Given: A set of "candidate" hi-DOFs to serve as the "representative" DOF for each lo-DOF on the reference element.
292 // Algorithm: For each element, we choose the lowest GID of the candidates for each DOF to generate the lo_elemToHiRepresentativeNode map
293
294 size_t numElem = hi_elemToNode.extent(0);
295 size_t lo_nperel = candidates.size();
296 Kokkos::resize(lo_elemToHiRepresentativeNode,numElem, lo_nperel);
297
298 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
299 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
300 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
301 for(size_t i=0; i<numElem; i++)
302 for(size_t j=0; j<lo_nperel; j++) {
303 if(candidates[j].size() == 1)
304 lo_elemToHiRepresentativeNode_host(i,j) = hi_elemToNode_host(i,candidates[j][0]);
305 else {
306 // First we get the GIDs for each candidate
307 std::vector<GO> GID(candidates[j].size());
308 for(size_t k=0; k<(size_t)candidates[j].size(); k++)
309 GID[k] = hi_columnMap->getGlobalElement(hi_elemToNode_host(i,candidates[j][k]));
310
311 // Find the one with smallest GID
312 size_t which = std::distance(GID.begin(),std::min_element(GID.begin(),GID.end()));
313
314 // Record this
315 lo_elemToHiRepresentativeNode_host(i,j) = hi_elemToNode_host(i,candidates[j][which]);
316 }
317 }
318 Kokkos::deep_copy(lo_elemToHiRepresentativeNode, lo_elemToHiRepresentativeNode_host);
319}
320
321/*********************************************************************************************************/
322// Inputs:
323// hi_elemToNode - FC<LO> containing the high order element-to-node map
324// hi_nodeIsOwned - std::vector<bool> of size hi's column map, which described hi node ownership
325// lo_elemToHiRepresentativeNode - FC<LO> of size (# elements, # lo dofs per element) listing the hi unknown chosen as the single representative for each lo unknown for counting purposes
326// Outputs:
327// lo_elemToNode - FC<LO> containing the low order element-to-node map.
328// lo_nodeIsOwned - std::vector<bool> of size lo's (future) column map, which described lo node ownership
329// hi_to_lo_map - std::vector<LO> of size equal to hi's column map, which contains the lo id each hi idea maps to (or invalid if it doesn't)
330// lo_numOwnedNodes- Number of lo owned nodes
331template <class LocalOrdinal, class LOFieldContainer>
332void BuildLoElemToNodeViaRepresentatives(const LOFieldContainer & hi_elemToNode,
333 const std::vector<bool> & hi_nodeIsOwned,
334 const LOFieldContainer & lo_elemToHiRepresentativeNode,
335 LOFieldContainer & lo_elemToNode,
336 std::vector<bool> & lo_nodeIsOwned,
337 std::vector<LocalOrdinal> & hi_to_lo_map,
338 int & lo_numOwnedNodes) {
339 typedef LocalOrdinal LO;
340 using Teuchos::RCP;
341 // printf("CMS:BuildLoElemToNodeViaRepresentatives: hi_elemToNode.rank() = %d hi_elemToNode.size() = %d\n",hi_elemToNode.rank(), hi_elemToNode.size());
342 size_t numElem = hi_elemToNode.extent(0);
343 size_t hi_numNodes = hi_nodeIsOwned.size();
344 size_t lo_nperel = lo_elemToHiRepresentativeNode.extent(1);
345 Kokkos::resize(lo_elemToNode,numElem, lo_nperel);
346
347 // Start by flagginc the representative nodes
348 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
349 Kokkos::deep_copy(lo_elemToHiRepresentativeNode_host, lo_elemToHiRepresentativeNode);
350 std::vector<bool> is_low_order(hi_numNodes,false);
351 for(size_t i=0; i<numElem; i++)
352 for(size_t j=0; j<lo_nperel; j++) {
353 LO id = lo_elemToHiRepresentativeNode_host(i,j);
354 is_low_order[id] = true; // This can overwrite and that is OK.
355 }
356
357 // Count the number of lo owned nodes, generating a local index for lo nodes
358 lo_numOwnedNodes=0;
359 size_t lo_numNodes=0;
360 hi_to_lo_map.resize(hi_numNodes,Teuchos::OrdinalTraits<LO>::invalid());
361
362 for(size_t i=0; i<hi_numNodes; i++)
363 if(is_low_order[i]) {
364 hi_to_lo_map[i] = lo_numNodes;
365 lo_numNodes++;
366 if(hi_nodeIsOwned[i]) lo_numOwnedNodes++;
367 }
368
369 // Flag the owned lo nodes
370 lo_nodeIsOwned.resize(lo_numNodes,false);
371 for(size_t i=0; i<hi_numNodes; i++) {
372 if(is_low_order[i] && hi_nodeIsOwned[i])
373 lo_nodeIsOwned[hi_to_lo_map[i]]=true;
374 }
375
376 // Translate lo_elemToNode to a lo local index
377 auto lo_elemToNode_host = Kokkos::create_mirror_view(lo_elemToNode);
378 for(size_t i=0; i<numElem; i++)
379 for(size_t j=0; j<lo_nperel; j++)
380 lo_elemToNode_host(i,j) = hi_to_lo_map[lo_elemToHiRepresentativeNode_host(i,j)];
381
382
383 // Check for the [E|T]petra column map ordering property, namely LIDs for owned nodes should all appear first.
384 // Since we're injecting from the higher-order mesh, it should be true, but we should add an error check & throw in case.
385 bool map_ordering_test_passed=true;
386 for(size_t i=0; i<lo_numNodes-1; i++)
387 if(!lo_nodeIsOwned[i] && lo_nodeIsOwned[i+1])
388 map_ordering_test_passed=false;
389
390 if(!map_ordering_test_passed)
391 throw std::runtime_error("MueLu::MueLuIntrepid::BuildLoElemToNodeViaRepresentatives failed map ordering test");
392 Kokkos::deep_copy(lo_elemToNode, lo_elemToNode_host);
393
394}
395
396
397/*********************************************************************************************************/
398// Inputs:
399// hi_elemToNode - FC<LO> containing the high order element-to-node map
400// hi_nodeIsOwned - std::vector<bool> of size hi's column map, which described hi node ownership
401// lo_node_in_hi - std::vector<size_t> of size lo dofs in the reference element, which describes the coindcident hi dots
402// hi_isDirichlet - ArrayView<int> of size of hi's column map, which has a 1 if the unknown is Dirichlet and a 0 if it isn't.
403// Outputs:
404// lo_elemToNode - FC<LO> containing the low order element-to-node map.
405// lo_nodeIsOwned - std::vector<bool> of size lo's (future) column map, which described lo node ownership
406// hi_to_lo_map - std::vector<LO> of size equal to hi's column map, which contains the lo id each hi idea maps to (or invalid if it doesn't)
407// lo_numOwnedNodes- Number of lo owned nodes
408template <class LocalOrdinal, class LOFieldContainer>
409void BuildLoElemToNode(const LOFieldContainer & hi_elemToNode,
410 const std::vector<bool> & hi_nodeIsOwned,
411 const std::vector<size_t> & lo_node_in_hi,
412 const Teuchos::ArrayRCP<const int> & hi_isDirichlet,
413 LOFieldContainer & lo_elemToNode,
414 std::vector<bool> & lo_nodeIsOwned,
415 std::vector<LocalOrdinal> & hi_to_lo_map,
416 int & lo_numOwnedNodes) {
417 typedef LocalOrdinal LO;
418 using Teuchos::RCP;
419 LocalOrdinal LOINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
420 // printf("CMS:BuildLoElemToNode: hi_elemToNode.rank() = %d hi_elemToNode.size() = %d\n",hi_elemToNode.rank(), hi_elemToNode.size());
421
422 size_t numElem = hi_elemToNode.extent(0);
423 size_t hi_numNodes = hi_nodeIsOwned.size();
424
425 size_t lo_nperel = lo_node_in_hi.size();
426 Kokkos::resize(lo_elemToNode,numElem, lo_nperel);
427
428 // Build lo_elemToNode (in the hi local index ordering) and flag owned ones
429 std::vector<bool> is_low_order(hi_numNodes,false);
430 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
431 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
432 auto lo_elemToNode_host = Kokkos::create_mirror_view(lo_elemToNode);
433 for(size_t i=0; i<numElem; i++)
434 for(size_t j=0; j<lo_nperel; j++) {
435 LO lid = hi_elemToNode_host(i,lo_node_in_hi[j]);
436
437 // Remove Dirichlet
438 if(hi_isDirichlet[lid])
439 lo_elemToNode_host(i,j) = LOINVALID;
440 else {
441 lo_elemToNode_host(i,j) = lid;
442 is_low_order[hi_elemToNode_host(i,lo_node_in_hi[j])] = true; // This can overwrite and that is OK.
443 }
444 }
445
446 // Count the number of lo owned nodes, generating a local index for lo nodes
447 lo_numOwnedNodes=0;
448 size_t lo_numNodes=0;
449 hi_to_lo_map.resize(hi_numNodes,Teuchos::OrdinalTraits<LO>::invalid());
450
451 for(size_t i=0; i<hi_numNodes; i++)
452 if(is_low_order[i]) {
453 hi_to_lo_map[i] = lo_numNodes;
454 lo_numNodes++;
455 if(hi_nodeIsOwned[i]) lo_numOwnedNodes++;
456 }
457
458 // Flag the owned lo nodes
459 lo_nodeIsOwned.resize(lo_numNodes,false);
460 for(size_t i=0; i<hi_numNodes; i++) {
461 if(is_low_order[i] && hi_nodeIsOwned[i])
462 lo_nodeIsOwned[hi_to_lo_map[i]]=true;
463 }
464
465 // Translate lo_elemToNode to a lo local index
466 for(size_t i=0; i<numElem; i++)
467 for(size_t j=0; j<lo_nperel; j++) {
468 if(lo_elemToNode_host(i,j) != LOINVALID)
469 lo_elemToNode_host(i,j) = hi_to_lo_map[lo_elemToNode_host(i,j)];
470 }
471 Kokkos::deep_copy(lo_elemToNode, lo_elemToNode_host);
472
473 // Check for the [E|T]petra column map ordering property, namely LIDs for owned nodes should all appear first.
474 // Since we're injecting from the higher-order mesh, it should be true, but we should add an error check & throw in case.
475 bool map_ordering_test_passed=true;
476 for(size_t i=0; i<lo_numNodes-1; i++)
477 if(!lo_nodeIsOwned[i] && lo_nodeIsOwned[i+1])
478 map_ordering_test_passed=false;
479
480 if(!map_ordering_test_passed)
481 throw std::runtime_error("MueLu::MueLuIntrepid::BuildLoElemToNode failed map ordering test");
482
483}
484
485/*********************************************************************************************************/
486// Generates the lo_columnMap
487// Input:
488// hi_importer - Importer from the hi matrix
489// hi_to_lo_map - std::vector<LO> of size equal to hi's column map, which contains the lo id each hi idea maps to (or invalid if it doesn't)
490// lo_DomainMap - Domain map for the lo matrix
491// lo_columnMapLength - Number of local columns in the lo column map
492// Output:
493// lo_columnMap - Column map of the lower order matrix
494 template <class LocalOrdinal, class GlobalOrdinal, class Node>
495 void GenerateColMapFromImport(const Xpetra::Import<LocalOrdinal,GlobalOrdinal, Node> & hi_importer,const std::vector<LocalOrdinal> &hi_to_lo_map,const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> & lo_domainMap, const size_t & lo_columnMapLength, RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > & lo_columnMap) {
496 typedef LocalOrdinal LO;
497 typedef GlobalOrdinal GO;
498 typedef Node NO;
499 typedef Xpetra::Map<LO,GO,NO> Map;
500 typedef Xpetra::Vector<GO,LO,GO,NO> GOVector;
501
502 GO go_invalid = Teuchos::OrdinalTraits<GO>::invalid();
503 LO lo_invalid = Teuchos::OrdinalTraits<LO>::invalid();
504
505 RCP<const Map> hi_domainMap = hi_importer.getSourceMap();
506 RCP<const Map> hi_columnMap = hi_importer.getTargetMap();
507 // Figure out the GIDs of my non-owned P1 nodes
508 // HOW: We can build a GOVector(domainMap) and fill the values with either invalid() or the P1 domainMap.GID() for that guy.
509 // Then we can use A's importer to get a GOVector(colMap) with that information.
510
511 // NOTE: This assumes rowMap==colMap and [E|T]petra ordering of all the locals first in the colMap
512 RCP<GOVector> dvec = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(hi_domainMap);
513 {
514 ArrayRCP<GO> dvec_data = dvec->getDataNonConst(0);
515 for(size_t i=0; i<hi_domainMap->getLocalNumElements(); i++) {
516 if(hi_to_lo_map[i]!=lo_invalid) dvec_data[i] = lo_domainMap.getGlobalElement(hi_to_lo_map[i]);
517 else dvec_data[i] = go_invalid;
518 }
519 }
520
521
522 RCP<GOVector> cvec = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(hi_columnMap,true);
523 cvec->doImport(*dvec,hi_importer,Xpetra::ADD);
524
525 // Generate the lo_columnMap
526 // HOW: We can use the local hi_to_lo_map from the GID's in cvec to generate the non-contiguous colmap ids.
527 Array<GO> lo_col_data(lo_columnMapLength);
528 {
529 ArrayRCP<GO> cvec_data = cvec->getDataNonConst(0);
530 for(size_t i=0,idx=0; i<hi_columnMap->getLocalNumElements(); i++) {
531 if(hi_to_lo_map[i]!=lo_invalid) {
532 lo_col_data[idx] = cvec_data[i];
533 idx++;
534 }
535 }
536 }
537
538 lo_columnMap = Xpetra::MapFactory<LO,GO,NO>::Build(lo_domainMap.lib(),Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),lo_col_data(),lo_domainMap.getIndexBase(),lo_domainMap.getComm());
539}
540
541/*********************************************************************************************************/
542// Generates a list of "representative candidate" hi dofs for each lo dof on the reference element. This is to be used in global numbering.
543// Input:
544// basis - The low order basis
545// ReferenceNodeLocations - FC<Scalar> of size (#hidofs, dim) Locations of higher order nodes on the reference element
546// threshold - tolerance for equivalance testing
547// Output:
548// representative_node_candidates - std::vector<std::vector<size_t> > of lists of "representative candidate" hi dofs for each lo dof
549template<class Basis, class SCFieldContainer>
550void GenerateRepresentativeBasisNodes(const Basis & basis, const SCFieldContainer & ReferenceNodeLocations, const double threshold,std::vector<std::vector<size_t> > & representative_node_candidates) {
551 typedef SCFieldContainer FC;
552 typedef typename FC::data_type SC;
553
554 // Evaluate the linear basis functions at the Pn nodes
555 size_t numFieldsHi = ReferenceNodeLocations.extent(0);
556 // size_t dim = ReferenceNodeLocations.extent(1);
557 size_t numFieldsLo = basis.getCardinality();
558
559 FC LoValues("LoValues",numFieldsLo,numFieldsHi);
560
561 basis.getValues(LoValues, ReferenceNodeLocations , Intrepid2::OPERATOR_VALUE);
562
563 Kokkos::fence(); // for kernel in getValues
564
565#if 0
566 printf("** LoValues[%d,%d] **\n",(int)numFieldsLo,(int)numFieldsHi);
567 for(size_t i=0; i<numFieldsLo; i++) {
568 for(size_t j=0; j<numFieldsHi; j++)
569 printf("%6.4e ",LoValues(i,j));
570 printf("\n");
571 }
572 printf("**************\n");fflush(stdout);
573#endif
574
575 representative_node_candidates.resize(numFieldsLo);
576 auto LoValues_host = Kokkos::create_mirror_view(LoValues);
577 Kokkos::deep_copy(LoValues_host, LoValues);
578 for(size_t i=0; i<numFieldsLo; i++) {
579 // 1st pass: find the max value
580 typename Teuchos::ScalarTraits<SC>::magnitudeType vmax = Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<SC>::magnitudeType>::zero();
581 for(size_t j=0; j<numFieldsHi; j++)
582 vmax = std::max(vmax,Teuchos::ScalarTraits<SC>::magnitude(LoValues_host(i,j)));
583
584 // 2nd pass: Find all values w/i threshhold of target
585 for(size_t j=0; j<numFieldsHi; j++) {
586 if(Teuchos::ScalarTraits<SC>::magnitude(vmax - LoValues_host(i,j)) < threshold*vmax)
587 representative_node_candidates[i].push_back(j);
588 }
589 }
590
591 // Sanity check
592 for(size_t i=0; i<numFieldsLo; i++)
593 if(!representative_node_candidates[i].size())
594 throw std::runtime_error("ERROR: GenerateRepresentativeBasisNodes: No candidates found!");
595
596
597}
598
599
600
601}//end MueLu::MueLuIntrepid namespace
602
603
604/*********************************************************************************************************/
605/*********************************************************************************************************/
606/*********************************************************************************************************/
607template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
609 const std::vector<bool> & hi_nodeIsOwned,
610 const SCFieldContainer & hi_DofCoords,
611 const std::vector<size_t> &lo_node_in_hi,
612 const Basis &lo_basis,
613 const std::vector<LocalOrdinal> & hi_to_lo_map,
614 const Teuchos::RCP<const Map> & lo_colMap,
615 const Teuchos::RCP<const Map> & lo_domainMap,
616 const Teuchos::RCP<const Map> & hi_map,
617 Teuchos::RCP<Matrix>& P) const{
618 typedef SCFieldContainer FC;
619 // Evaluate the linear basis functions at the Pn nodes
620 size_t numFieldsHi = hi_elemToNode.extent(1);
621 size_t numFieldsLo = lo_basis.getCardinality();
622 LocalOrdinal LOINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
623 FC LoValues_at_HiDofs("LoValues_at_HiDofs",numFieldsLo,numFieldsHi);
624 lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
625 auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
626 Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
627 Kokkos::fence(); // for kernel in getValues
628
629 typedef typename Teuchos::ScalarTraits<SC>::halfPrecision SClo;
630 typedef typename Teuchos::ScalarTraits<SClo>::magnitudeType MT;
631 MT effective_zero = Teuchos::ScalarTraits<MT>::eps();
632
633 // Allocate P
634 P = rcp(new CrsMatrixWrap(hi_map,lo_colMap,numFieldsHi)); //FIXLATER: Need faster fill
635 RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
636
637 // Slow-ish fill
638 size_t Nelem=hi_elemToNode.extent(0);
639 std::vector<bool> touched(hi_map->getLocalNumElements(),false);
640 Teuchos::Array<GO> col_gid(1);
641 Teuchos::Array<SC> val(1);
642 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
643 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
644 for(size_t i=0; i<Nelem; i++) {
645 for(size_t j=0; j<numFieldsHi; j++) {
646 LO row_lid = hi_elemToNode_host(i,j);
647 GO row_gid = hi_map->getGlobalElement(row_lid);
648 if(hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
649 for(size_t k=0; k<numFieldsLo; k++) {
650 // Get the local id in P1's column map
651 LO col_lid = hi_to_lo_map[hi_elemToNode_host(i,lo_node_in_hi[k])];
652 if(col_lid==LOINVALID) continue;
653
654 col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
655 val[0] = LoValues_at_HiDofs_host(k,j);
656
657 // Skip near-zeros
658 if(Teuchos::ScalarTraits<SC>::magnitude(val[0]) >= effective_zero)
659 P->insertGlobalValues(row_gid,col_gid(),val());
660 }
661 touched[row_lid]=true;
662 }
663 }
664 }
665 P->fillComplete(lo_domainMap,hi_map);
666}
667
668/*********************************************************************************************************/
669template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
671 const std::vector<bool> & hi_nodeIsOwned,
672 const SCFieldContainer & hi_DofCoords,
673 const LOFieldContainer & lo_elemToHiRepresentativeNode,
674 const Basis &lo_basis,
675 const std::vector<LocalOrdinal> & hi_to_lo_map,
676 const Teuchos::RCP<const Map> & lo_colMap,
677 const Teuchos::RCP<const Map> & lo_domainMap,
678 const Teuchos::RCP<const Map> & hi_map,
679 Teuchos::RCP<Matrix>& P) const{
680 typedef SCFieldContainer FC;
681 // Evaluate the linear basis functions at the Pn nodes
682 size_t numFieldsHi = hi_elemToNode.extent(1);
683 size_t numFieldsLo = lo_basis.getCardinality();
684 FC LoValues_at_HiDofs("LoValues_at_HiDofs",numFieldsLo,numFieldsHi);
685 lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
686 auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
687 auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
688 auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
689 Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
690 Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
691 Kokkos::deep_copy(lo_elemToHiRepresentativeNode_host, lo_elemToHiRepresentativeNode);
692 Kokkos::fence(); // for kernel in getValues
693
694 typedef typename Teuchos::ScalarTraits<SC>::halfPrecision SClo;
695 typedef typename Teuchos::ScalarTraits<SClo>::magnitudeType MT;
696 MT effective_zero = Teuchos::ScalarTraits<MT>::eps();
697
698 // Allocate P
699 P = rcp(new CrsMatrixWrap(hi_map,lo_colMap,numFieldsHi)); //FIXLATER: Need faster fill
700 RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
701
702 // Slow-ish fill
703 size_t Nelem=hi_elemToNode.extent(0);
704 std::vector<bool> touched(hi_map->getLocalNumElements(),false);
705 Teuchos::Array<GO> col_gid(1);
706 Teuchos::Array<SC> val(1);
707 for(size_t i=0; i<Nelem; i++) {
708 for(size_t j=0; j<numFieldsHi; j++) {
709 LO row_lid = hi_elemToNode_host(i,j);
710 GO row_gid = hi_map->getGlobalElement(row_lid);
711 if(hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
712 for(size_t k=0; k<numFieldsLo; k++) {
713 // Get the local id in P1's column map
714 LO col_lid = hi_to_lo_map[lo_elemToHiRepresentativeNode_host(i,k)];
715 col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
716 val[0] = LoValues_at_HiDofs_host(k,j);
717
718 // Skip near-zeros
719 if(Teuchos::ScalarTraits<SC>::magnitude(val[0]) >= effective_zero)
720 P->insertGlobalValues(row_gid,col_gid(),val());
721 }
722 touched[row_lid]=true;
723 }
724 }
725 }
726 P->fillComplete(lo_domainMap,hi_map);
727}
728
729/*********************************************************************************************************/
730 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
732 RCP<ParameterList> validParamList = rcp(new ParameterList());
733
734#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
735 SET_VALID_ENTRY("pcoarsen: hi basis");
736 SET_VALID_ENTRY("pcoarsen: lo basis");
737#undef SET_VALID_ENTRY
738
739 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A used during the prolongator smoothing process");
740
741 validParamList->set< RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory of the nullspace");
742 validParamList->set< RCP<const FactoryBase> >("pcoarsen: element to node map", Teuchos::null, "Generating factory of the element to node map");
743 return validParamList;
744 }
745
746/*********************************************************************************************************/
747 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
749 Input(fineLevel, "A");
750 Input(fineLevel, "pcoarsen: element to node map");
751 Input(fineLevel, "Nullspace");
752 }
753
754/*********************************************************************************************************/
755 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
757 return BuildP(fineLevel, coarseLevel);
758 }
759
760/*********************************************************************************************************/
761 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
763 FactoryMonitor m(*this, "P Coarsening", coarseLevel);
764 std::string levelIDs = toString(coarseLevel.GetLevelID());
765 const std::string prefix = "MueLu::IntrepidPCoarsenFactory(" + levelIDs + "): ";
766
767 // NOTE: This is hardwired to double on purpose. See the note below.
768 typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCi;
769 typedef Kokkos::DynRankView<double,typename Node::device_type> FC;
770
771 // Level Get
772 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
773 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> >(fineLevel, "Nullspace");
774 Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Acrs = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*A);
775
776
777 if (restrictionMode_) {
778 SubFactoryMonitor m2(*this, "Transpose A", coarseLevel);
779 A = Utilities::Transpose(*A, true); // build transpose of A explicitly
780 }
781
782 // Find the Dirichlet rows in A
783 std::vector<LocalOrdinal> A_dirichletRows;
784 Utilities::FindDirichletRows(A,A_dirichletRows);
785
786 // Build final prolongator
787 RCP<Matrix> finalP;
788
789 // Reuse pattern if available
790 RCP<ParameterList> APparams = rcp(new ParameterList);
791 if (coarseLevel.IsAvailable("AP reuse data", this)) {
792 GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous AP data" << std::endl;
793
794 APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
795
796 if (APparams->isParameter("graph"))
797 finalP = APparams->get< RCP<Matrix> >("graph");
798 }
799 const ParameterList& pL = GetParameterList();
800
801 /*******************/
802 // FIXME LATER: Allow these to be manually specified instead of Intrepid
803 // Get the Intrepid bases
804 // NOTE: To make sure Stokhos works we only instantiate these guys with double. There's a lot
805 // of stuff in the guts of Intrepid2 that doesn't play well with Stokhos as of yet.
806 int lo_degree, hi_degree;
807 RCP<Basis> hi_basis = MueLuIntrepid::BasisFactory<double,typename Node::device_type::execution_space>(pL.get<std::string>("pcoarsen: hi basis"),hi_degree);
808 RCP<Basis> lo_basis = MueLuIntrepid::BasisFactory<double,typename Node::device_type::execution_space>(pL.get<std::string>("pcoarsen: lo basis"),lo_degree);
809
810 // Useful Output
811 GetOStream(Statistics1) << "P-Coarsening from basis "<<pL.get<std::string>("pcoarsen: hi basis")<<" to "<<pL.get<std::string>("pcoarsen: lo basis") <<std::endl;
812
813 /*******************/
814 // Get the higher-order element-to-node map
815 const Teuchos::RCP<FCi> Pn_elemToNode = Get<Teuchos::RCP<FCi> >(fineLevel,"pcoarsen: element to node map");
816
817 // Calculate node ownership (the quick and dirty way)
818 // NOTE: This exploits two things:
819 // 1) domainMap == rowMap
820 // 2) Standard [e|t]petra ordering (namely the local unknowns are always numbered first).
821 // This routine does not work in general.
822 RCP<const Map> rowMap = A->getRowMap();
823 RCP<const Map> colMap = Acrs.getColMap();
824 RCP<const Map> domainMap = A->getDomainMap();
825 int NumProc = rowMap->getComm()->getSize();
826 assert(rowMap->isSameAs(*domainMap));
827 std::vector<bool> Pn_nodeIsOwned(colMap->getLocalNumElements(),false);
828 LO num_owned_rows = 0;
829 for(size_t i=0; i<rowMap->getLocalNumElements(); i++) {
830 if(rowMap->getGlobalElement(i) == colMap->getGlobalElement(i)) {
831 Pn_nodeIsOwned[i] = true;
832 num_owned_rows++;
833 }
834 }
835
836 // Used in all cases
837 FC hi_DofCoords;
838 Teuchos::RCP<FCi> P1_elemToNode = rcp(new FCi());
839
840 std::vector<bool> P1_nodeIsOwned;
841 int P1_numOwnedNodes;
842 std::vector<LO> hi_to_lo_map;
843
844 // Degree-1 variables
845 std::vector<size_t> lo_node_in_hi;
846
847 // Degree-n variables
848 FCi lo_elemToHiRepresentativeNode;
849
850 // Get Dirichlet unknown information
851 RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> > hi_isDirichletRow, hi_isDirichletCol;
852 Utilities::FindDirichletRowsAndPropagateToCols(A,hi_isDirichletRow, hi_isDirichletCol);
853
854#if 0
855 printf("[%d] isDirichletRow = ",A->getRowMap()->getComm()->getRank());
856 for(size_t i=0;i<hi_isDirichletRow->getMap()->getLocalNumElements(); i++)
857 printf("%d ",hi_isDirichletRow->getData(0)[i]);
858 printf("\n");
859 printf("[%d] isDirichletCol = ",A->getRowMap()->getComm()->getRank());
860 for(size_t i=0;i<hi_isDirichletCol->getMap()->getLocalNumElements(); i++)
861 printf("%d ",hi_isDirichletCol->getData(0)[i]);
862 printf("\n");
863 fflush(stdout);
864#endif
865
866 /*******************/
867 if(lo_degree == 1) {
868 // Get reference coordinates and the lo-to-hi injection list for the reference element
869 MueLuIntrepid::IntrepidGetP1NodeInHi(hi_basis,lo_node_in_hi,hi_DofCoords);
870
871 // Generate lower-order element-to-node map
872 MueLuIntrepid::BuildLoElemToNode(*Pn_elemToNode,Pn_nodeIsOwned,lo_node_in_hi,hi_isDirichletCol->getData(0),*P1_elemToNode,P1_nodeIsOwned,hi_to_lo_map,P1_numOwnedNodes);
873 assert(hi_to_lo_map.size() == colMap->getLocalNumElements());
874 }
875 else {
876 // Get lo-order candidates
877 double threshold = 1e-10;
878 std::vector<std::vector<size_t> > candidates;
879 Kokkos::resize(hi_DofCoords,hi_basis->getCardinality(),hi_basis->getBaseCellTopology().getDimension());
880 hi_basis->getDofCoords(hi_DofCoords);
881
882 MueLu::MueLuIntrepid::GenerateRepresentativeBasisNodes<Basis,FC>(*lo_basis,hi_DofCoords,threshold,candidates);
883
884 // Generate the representative nodes
885 MueLu::MueLuIntrepid::GenerateLoNodeInHiViaGIDs(candidates,*Pn_elemToNode,colMap,lo_elemToHiRepresentativeNode);
886 MueLu::MueLuIntrepid::BuildLoElemToNodeViaRepresentatives(*Pn_elemToNode,Pn_nodeIsOwned,lo_elemToHiRepresentativeNode,*P1_elemToNode,P1_nodeIsOwned,hi_to_lo_map,P1_numOwnedNodes);
887 }
888 MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(coarseLevel,"pcoarsen: element to node map",P1_elemToNode);
889
890 /*******************/
891 // Generate the P1_domainMap
892 // HOW: Since we know how many each proc has, we can use the non-uniform contiguous map constructor to do the work for us
893 RCP<const Map> P1_domainMap = MapFactory::Build(rowMap->lib(),Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),P1_numOwnedNodes,rowMap->getIndexBase(),rowMap->getComm());
894 MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(coarseLevel,"CoarseMap",P1_domainMap);
895
896 // Generate the P1_columnMap
897 RCP<const Map> P1_colMap;
898 if(NumProc==1) P1_colMap = P1_domainMap;
899 else MueLuIntrepid::GenerateColMapFromImport<LO,GO,NO>(*Acrs.getCrsGraph()->getImporter(),hi_to_lo_map,*P1_domainMap,P1_nodeIsOwned.size(),P1_colMap);
900
901 /*******************/
902 // Generate the coarsening
903 if(lo_degree == 1)
904 GenerateLinearCoarsening_pn_kirby_to_p1(*Pn_elemToNode,Pn_nodeIsOwned,hi_DofCoords,lo_node_in_hi,*lo_basis,hi_to_lo_map,P1_colMap,P1_domainMap,A->getRowMap(),finalP);
905 else
906 GenerateLinearCoarsening_pn_kirby_to_pm(*Pn_elemToNode,Pn_nodeIsOwned,hi_DofCoords,lo_elemToHiRepresentativeNode,*lo_basis,hi_to_lo_map,P1_colMap,P1_domainMap,A->getRowMap(),finalP);
907
908 /*******************/
909 // Zero out the Dirichlet rows in P
910 Utilities::ZeroDirichletRows(finalP,A_dirichletRows);
911
912 /*******************/
913 // Build the nullspace
914 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P1_domainMap, fineNullspace->getNumVectors());
915 finalP->apply(*fineNullspace,*coarseNullspace,Teuchos::TRANS);
916 Set(coarseLevel, "Nullspace", coarseNullspace);
917
918 // Level Set
919 if (!restrictionMode_) {
920 // The factory is in prolongation mode
921 Set(coarseLevel, "P", finalP);
922
923 APparams->set("graph", finalP);
924 MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(coarseLevel,"AP reuse data",APparams);
925
926 if (IsPrint(Statistics1)) {
927 RCP<ParameterList> params = rcp(new ParameterList());
928 params->set("printLoadBalancingInfo", true);
929 params->set("printCommInfo", true);
930 GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*finalP, "P", params);
931 }
932 } else {
933 // The factory is in restriction mode
934 RCP<Matrix> R;
935 {
936 SubFactoryMonitor m2(*this, "Transpose P", coarseLevel);
937 R = Utilities::Transpose(*finalP, true);
938 }
939
940 Set(coarseLevel, "R", R);
941
942 if (IsPrint(Statistics2)) {
943 RCP<ParameterList> params = rcp(new ParameterList());
944 params->set("printLoadBalancingInfo", true);
945 params->set("printCommInfo", true);
946 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*R, "R", params);
947 }
948 }
949
950 } //Build()
951
952} //namespace MueLu
953
954#endif // MUELU_IPCFACTORY_DEF_HPP
955
#define SET_VALID_ENTRY(name)
#define MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(level, ename, entry)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
Timer to be used in factories. Similar to Monitor but with additional timers.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Kokkos::DynRankView< LocalOrdinal, typename Node::device_type > LOFieldContainer
void GenerateLinearCoarsening_pn_kirby_to_p1(const LOFieldContainer &hi_elemToNode, const std::vector< bool > &hi_nodeIsOwned, const SCFieldContainer &hi_DofCoords, const std::vector< size_t > &lo_node_in_hi, const Basis &lo_Basis, const std::vector< LocalOrdinal > &hi_to_lo_map, const Teuchos::RCP< const Map > &lo_colMap, const Teuchos::RCP< const Map > &lo_domainMap, const Teuchos::RCP< const Map > &hi_map, Teuchos::RCP< Matrix > &P) const
Kokkos::DynRankView< double, typename Node::device_type > SCFieldContainer
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Intrepid2::Basis< typename Node::device_type::execution_space, double, double > Basis
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void GenerateLinearCoarsening_pn_kirby_to_pm(const LOFieldContainer &hi_elemToNode, const std::vector< bool > &hi_nodeIsOwned, const SCFieldContainer &hi_DofCoords, const LOFieldContainer &lo_elemToHiRepresentativeNode, const Basis &lo_basis, const std::vector< LocalOrdinal > &hi_to_lo_map, const Teuchos::RCP< const Map > &lo_colMap, const Teuchos::RCP< const Map > &lo_domainMap, const Teuchos::RCP< const Map > &hi_map, Teuchos::RCP< Matrix > &P) const
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node > > &isDirichletRow, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node > > &isDirichletCol)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
void IntrepidGetP1NodeInHi(const Teuchos::RCP< Intrepid2::Basis< typename KokkosDeviceType::execution_space, Scalar, Scalar > > &hi_basis, std::vector< size_t > &lo_node_in_hi, Kokkos::DynRankView< Scalar, KokkosDeviceType > &hi_DofCoords)
void FindGeometricSeedOrdinals(Teuchos::RCP< Basis > basis, const LOFieldContainer &elementToNodeMap, std::vector< std::vector< LocalOrdinal > > &seeds, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &columnMap)
void GenerateRepresentativeBasisNodes(const Basis &basis, const SCFieldContainer &ReferenceNodeLocations, const double threshold, std::vector< std::vector< size_t > > &representative_node_candidates)
void GenerateColMapFromImport(const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &hi_importer, const std::vector< LocalOrdinal > &hi_to_lo_map, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &lo_domainMap, const size_t &lo_columnMapLength, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &lo_columnMap)
void BuildLoElemToNodeViaRepresentatives(const LOFieldContainer &hi_elemToNode, const std::vector< bool > &hi_nodeIsOwned, const LOFieldContainer &lo_elemToHiRepresentativeNode, LOFieldContainer &lo_elemToNode, std::vector< bool > &lo_nodeIsOwned, std::vector< LocalOrdinal > &hi_to_lo_map, int &lo_numOwnedNodes)
std::string tolower(const std::string &str)
void BuildLoElemToNode(const LOFieldContainer &hi_elemToNode, const std::vector< bool > &hi_nodeIsOwned, const std::vector< size_t > &lo_node_in_hi, const Teuchos::ArrayRCP< const int > &hi_isDirichlet, LOFieldContainer &lo_elemToNode, std::vector< bool > &lo_nodeIsOwned, std::vector< LocalOrdinal > &hi_to_lo_map, int &lo_numOwnedNodes)
void GenerateLoNodeInHiViaGIDs(const std::vector< std::vector< size_t > > &candidates, const LOFieldContainer &hi_elemToNode, RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &hi_columnMap, LOFieldContainer &lo_elemToHiRepresentativeNode)
Teuchos::RCP< Intrepid2::Basis< KokkosExecutionSpace, Scalar, Scalar > > BasisFactory(const std::string &name, int &degree)
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Statistics1
Print more statistics.
@ Runtime0
One-liner description of what is happening.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.