MueLu  Version of the Day
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 
100 namespace MueLu {
101 
102 
103 /*********************************************************************************************************/
104 namespace MueLuIntrepid {
105 inline 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
196 template<class Scalar,class KokkosExecutionSpace>
197 Teuchos::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
250 template<class Scalar,class KokkosDeviceType>
251 void 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
285 template<class LocalOrdinal, class GlobalOrdinal, class Node, class LOFieldContainer>
286 void 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
331 template <class LocalOrdinal, class LOFieldContainer>
332 void 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
408 template <class LocalOrdinal, class LOFieldContainer>
409 void 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  ArrayRCP<GO> dvec_data = dvec->getDataNonConst(0);
514  for(size_t i=0; i<hi_domainMap->getNodeNumElements(); i++) {
515  if(hi_to_lo_map[i]!=lo_invalid) dvec_data[i] = lo_domainMap.getGlobalElement(hi_to_lo_map[i]);
516  else dvec_data[i] = go_invalid;
517  }
518 
519 
520  RCP<GOVector> cvec = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(hi_columnMap,true);
521  cvec->doImport(*dvec,hi_importer,Xpetra::ADD);
522 
523  // Generate the lo_columnMap
524  // HOW: We can use the local hi_to_lo_map from the GID's in cvec to generate the non-contiguous colmap ids.
525  Array<GO> lo_col_data(lo_columnMapLength);
526  ArrayRCP<GO> cvec_data = cvec->getDataNonConst(0);
527  for(size_t i=0,idx=0; i<hi_columnMap->getNodeNumElements(); i++) {
528  if(hi_to_lo_map[i]!=lo_invalid) {
529  lo_col_data[idx] = cvec_data[i];
530  idx++;
531  }
532  }
533 
534  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());
535 }
536 
537 /*********************************************************************************************************/
538 // Generates a list of "representative candidate" hi dofs for each lo dof on the reference element. This is to be used in global numbering.
539 // Input:
540 // basis - The low order basis
541 // ReferenceNodeLocations - FC<Scalar> of size (#hidofs, dim) Locations of higher order nodes on the reference element
542 // threshold - tolerance for equivalance testing
543 // Output:
544 // representative_node_candidates - std::vector<std::vector<size_t> > of lists of "representative candidate" hi dofs for each lo dof
545 template<class Basis, class SCFieldContainer>
546 void GenerateRepresentativeBasisNodes(const Basis & basis, const SCFieldContainer & ReferenceNodeLocations, const double threshold,std::vector<std::vector<size_t> > & representative_node_candidates) {
547  typedef SCFieldContainer FC;
548  typedef typename FC::data_type SC;
549 
550  // Evaluate the linear basis functions at the Pn nodes
551  size_t numFieldsHi = ReferenceNodeLocations.extent(0);
552  // size_t dim = ReferenceNodeLocations.extent(1);
553  size_t numFieldsLo = basis.getCardinality();
554 
555  FC LoValues("LoValues",numFieldsLo,numFieldsHi);
556 
557  basis.getValues(LoValues, ReferenceNodeLocations , Intrepid2::OPERATOR_VALUE);
558 
559  Kokkos::fence(); // for kernel in getValues
560 
561 #if 0
562  printf("** LoValues[%d,%d] **\n",(int)numFieldsLo,(int)numFieldsHi);
563  for(size_t i=0; i<numFieldsLo; i++) {
564  for(size_t j=0; j<numFieldsHi; j++)
565  printf("%6.4e ",LoValues(i,j));
566  printf("\n");
567  }
568  printf("**************\n");fflush(stdout);
569 #endif
570 
571  representative_node_candidates.resize(numFieldsLo);
572  auto LoValues_host = Kokkos::create_mirror_view(LoValues);
573  Kokkos::deep_copy(LoValues_host, LoValues);
574  for(size_t i=0; i<numFieldsLo; i++) {
575  // 1st pass: find the max value
576  typename Teuchos::ScalarTraits<SC>::magnitudeType vmax = Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<SC>::magnitudeType>::zero();
577  for(size_t j=0; j<numFieldsHi; j++)
578  vmax = std::max(vmax,Teuchos::ScalarTraits<SC>::magnitude(LoValues_host(i,j)));
579 
580  // 2nd pass: Find all values w/i threshhold of target
581  for(size_t j=0; j<numFieldsHi; j++) {
582  if(Teuchos::ScalarTraits<SC>::magnitude(vmax - LoValues_host(i,j)) < threshold*vmax)
583  representative_node_candidates[i].push_back(j);
584  }
585  }
586 
587  // Sanity check
588  for(size_t i=0; i<numFieldsLo; i++)
589  if(!representative_node_candidates[i].size())
590  throw std::runtime_error("ERROR: GenerateRepresentativeBasisNodes: No candidates found!");
591 
592 
593 }
594 
595 
596 
597 }//end MueLu::MueLuIntrepid namespace
598 
599 
600 /*********************************************************************************************************/
601 /*********************************************************************************************************/
602 /*********************************************************************************************************/
603 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
605  const std::vector<bool> & hi_nodeIsOwned,
606  const SCFieldContainer & hi_DofCoords,
607  const std::vector<size_t> &lo_node_in_hi,
608  const Basis &lo_basis,
609  const std::vector<LocalOrdinal> & hi_to_lo_map,
610  const Teuchos::RCP<const Map> & lo_colMap,
611  const Teuchos::RCP<const Map> & lo_domainMap,
612  const Teuchos::RCP<const Map> & hi_map,
613  Teuchos::RCP<Matrix>& P) const{
614  typedef SCFieldContainer FC;
615  // Evaluate the linear basis functions at the Pn nodes
616  size_t numFieldsHi = hi_elemToNode.extent(1);
617  size_t numFieldsLo = lo_basis.getCardinality();
618  LocalOrdinal LOINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
619  FC LoValues_at_HiDofs("LoValues_at_HiDofs",numFieldsLo,numFieldsHi);
620  lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
621  auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
622  Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
623  Kokkos::fence(); // for kernel in getValues
624 
625  typedef typename Teuchos::ScalarTraits<SC>::halfPrecision SClo;
626  typedef typename Teuchos::ScalarTraits<SClo>::magnitudeType MT;
627  MT effective_zero = Teuchos::ScalarTraits<MT>::eps();
628 
629  // Allocate P
630  P = rcp(new CrsMatrixWrap(hi_map,lo_colMap,numFieldsHi)); //FIXLATER: Need faster fill
631  RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
632 
633  // Slow-ish fill
634  size_t Nelem=hi_elemToNode.extent(0);
635  std::vector<bool> touched(hi_map->getNodeNumElements(),false);
636  Teuchos::Array<GO> col_gid(1);
637  Teuchos::Array<SC> val(1);
638  auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
639  Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
640  for(size_t i=0; i<Nelem; i++) {
641  for(size_t j=0; j<numFieldsHi; j++) {
642  LO row_lid = hi_elemToNode_host(i,j);
643  GO row_gid = hi_map->getGlobalElement(row_lid);
644  if(hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
645  for(size_t k=0; k<numFieldsLo; k++) {
646  // Get the local id in P1's column map
647  LO col_lid = hi_to_lo_map[hi_elemToNode_host(i,lo_node_in_hi[k])];
648  if(col_lid==LOINVALID) continue;
649 
650  col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
651  val[0] = LoValues_at_HiDofs_host(k,j);
652 
653  // Skip near-zeros
654  if(Teuchos::ScalarTraits<SC>::magnitude(val[0]) >= effective_zero)
655  P->insertGlobalValues(row_gid,col_gid(),val());
656  }
657  touched[row_lid]=true;
658  }
659  }
660  }
661  P->fillComplete(lo_domainMap,hi_map);
662 }
663 
664 /*********************************************************************************************************/
665 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
667  const std::vector<bool> & hi_nodeIsOwned,
668  const SCFieldContainer & hi_DofCoords,
669  const LOFieldContainer & lo_elemToHiRepresentativeNode,
670  const Basis &lo_basis,
671  const std::vector<LocalOrdinal> & hi_to_lo_map,
672  const Teuchos::RCP<const Map> & lo_colMap,
673  const Teuchos::RCP<const Map> & lo_domainMap,
674  const Teuchos::RCP<const Map> & hi_map,
675  Teuchos::RCP<Matrix>& P) const{
676  typedef SCFieldContainer FC;
677  // Evaluate the linear basis functions at the Pn nodes
678  size_t numFieldsHi = hi_elemToNode.extent(1);
679  size_t numFieldsLo = lo_basis.getCardinality();
680  FC LoValues_at_HiDofs("LoValues_at_HiDofs",numFieldsLo,numFieldsHi);
681  lo_basis.getValues(LoValues_at_HiDofs, hi_DofCoords, Intrepid2::OPERATOR_VALUE);
682  auto LoValues_at_HiDofs_host = Kokkos::create_mirror_view(LoValues_at_HiDofs);
683  auto hi_elemToNode_host = Kokkos::create_mirror_view(hi_elemToNode);
684  auto lo_elemToHiRepresentativeNode_host = Kokkos::create_mirror_view(lo_elemToHiRepresentativeNode);
685  Kokkos::deep_copy(LoValues_at_HiDofs_host, LoValues_at_HiDofs);
686  Kokkos::deep_copy(hi_elemToNode_host, hi_elemToNode);
687  Kokkos::deep_copy(lo_elemToHiRepresentativeNode_host, lo_elemToHiRepresentativeNode);
688  Kokkos::fence(); // for kernel in getValues
689 
690  typedef typename Teuchos::ScalarTraits<SC>::halfPrecision SClo;
691  typedef typename Teuchos::ScalarTraits<SClo>::magnitudeType MT;
692  MT effective_zero = Teuchos::ScalarTraits<MT>::eps();
693 
694  // Allocate P
695  P = rcp(new CrsMatrixWrap(hi_map,lo_colMap,numFieldsHi)); //FIXLATER: Need faster fill
696  RCP<CrsMatrix> Pcrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
697 
698  // Slow-ish fill
699  size_t Nelem=hi_elemToNode.extent(0);
700  std::vector<bool> touched(hi_map->getNodeNumElements(),false);
701  Teuchos::Array<GO> col_gid(1);
702  Teuchos::Array<SC> val(1);
703  for(size_t i=0; i<Nelem; i++) {
704  for(size_t j=0; j<numFieldsHi; j++) {
705  LO row_lid = hi_elemToNode_host(i,j);
706  GO row_gid = hi_map->getGlobalElement(row_lid);
707  if(hi_nodeIsOwned[row_lid] && !touched[row_lid]) {
708  for(size_t k=0; k<numFieldsLo; k++) {
709  // Get the local id in P1's column map
710  LO col_lid = hi_to_lo_map[lo_elemToHiRepresentativeNode_host(i,k)];
711  col_gid[0] = {lo_colMap->getGlobalElement(col_lid)};
712  val[0] = LoValues_at_HiDofs_host(k,j);
713 
714  // Skip near-zeros
715  if(Teuchos::ScalarTraits<SC>::magnitude(val[0]) >= effective_zero)
716  P->insertGlobalValues(row_gid,col_gid(),val());
717  }
718  touched[row_lid]=true;
719  }
720  }
721  }
722  P->fillComplete(lo_domainMap,hi_map);
723 }
724 
725 /*********************************************************************************************************/
726  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
728  RCP<ParameterList> validParamList = rcp(new ParameterList());
729 
730 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
731  SET_VALID_ENTRY("pcoarsen: hi basis");
732  SET_VALID_ENTRY("pcoarsen: lo basis");
733 #undef SET_VALID_ENTRY
734 
735  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A used during the prolongator smoothing process");
736 
737  validParamList->set< RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory of the nullspace");
738  validParamList->set< RCP<const FactoryBase> >("pcoarsen: element to node map", Teuchos::null, "Generating factory of the element to node map");
739  return validParamList;
740  }
741 
742 /*********************************************************************************************************/
743  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
745  Input(fineLevel, "A");
746  Input(fineLevel, "pcoarsen: element to node map");
747  Input(fineLevel, "Nullspace");
748  }
749 
750 /*********************************************************************************************************/
751  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
753  return BuildP(fineLevel, coarseLevel);
754  }
755 
756 /*********************************************************************************************************/
757  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
759  FactoryMonitor m(*this, "P Coarsening", coarseLevel);
760  std::string levelIDs = toString(coarseLevel.GetLevelID());
761  const std::string prefix = "MueLu::IntrepidPCoarsenFactory(" + levelIDs + "): ";
762 
763  // NOTE: This is hardwired to double on purpose. See the note below.
764  typedef Kokkos::DynRankView<LocalOrdinal,typename Node::device_type> FCi;
765  typedef Kokkos::DynRankView<double,typename Node::device_type> FC;
766 
767  // Level Get
768  RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
769  RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> >(fineLevel, "Nullspace");
770  Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Acrs = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*A);
771 
772 
773  if (restrictionMode_) {
774  SubFactoryMonitor m2(*this, "Transpose A", coarseLevel);
775  A = Utilities::Transpose(*A, true); // build transpose of A explicitly
776  }
777 
778  // Find the Dirichlet rows in A
779  std::vector<LocalOrdinal> A_dirichletRows;
780  Utilities::FindDirichletRows(A,A_dirichletRows);
781 
782  // Build final prolongator
783  RCP<Matrix> finalP;
784 
785  // Reuse pattern if available
786  RCP<ParameterList> APparams = rcp(new ParameterList);
787  if (coarseLevel.IsAvailable("AP reuse data", this)) {
788  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous AP data" << std::endl;
789 
790  APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
791 
792  if (APparams->isParameter("graph"))
793  finalP = APparams->get< RCP<Matrix> >("graph");
794  }
795  const ParameterList& pL = GetParameterList();
796 
797  /*******************/
798  // FIXME LATER: Allow these to be manually specified instead of Intrepid
799  // Get the Intrepid bases
800  // NOTE: To make sure Stokhos works we only instantiate these guys with double. There's a lot
801  // of stuff in the guts of Intrepid2 that doesn't play well with Stokhos as of yet.
802  int lo_degree, hi_degree;
803  RCP<Basis> hi_basis = MueLuIntrepid::BasisFactory<double,typename Node::device_type::execution_space>(pL.get<std::string>("pcoarsen: hi basis"),hi_degree);
804  RCP<Basis> lo_basis = MueLuIntrepid::BasisFactory<double,typename Node::device_type::execution_space>(pL.get<std::string>("pcoarsen: lo basis"),lo_degree);
805 
806  // Useful Output
807  GetOStream(Statistics1) << "P-Coarsening from basis "<<pL.get<std::string>("pcoarsen: hi basis")<<" to "<<pL.get<std::string>("pcoarsen: lo basis") <<std::endl;
808 
809  /*******************/
810  // Get the higher-order element-to-node map
811  const Teuchos::RCP<FCi> Pn_elemToNode = Get<Teuchos::RCP<FCi> >(fineLevel,"pcoarsen: element to node map");
812 
813  // Calculate node ownership (the quick and dirty way)
814  // NOTE: This exploits two things:
815  // 1) domainMap == rowMap
816  // 2) Standard [e|t]petra ordering (namely the local unknowns are always numbered first).
817  // This routine does not work in general.
818  RCP<const Map> rowMap = A->getRowMap();
819  RCP<const Map> colMap = Acrs.getColMap();
820  RCP<const Map> domainMap = A->getDomainMap();
821  int NumProc = rowMap->getComm()->getSize();
822  assert(rowMap->isSameAs(*domainMap));
823  std::vector<bool> Pn_nodeIsOwned(colMap->getNodeNumElements(),false);
824  LO num_owned_rows = 0;
825  for(size_t i=0; i<rowMap->getNodeNumElements(); i++) {
826  if(rowMap->getGlobalElement(i) == colMap->getGlobalElement(i)) {
827  Pn_nodeIsOwned[i] = true;
828  num_owned_rows++;
829  }
830  }
831 
832  // Used in all cases
833  FC hi_DofCoords;
834  Teuchos::RCP<FCi> P1_elemToNode = rcp(new FCi());
835 
836  std::vector<bool> P1_nodeIsOwned;
837  int P1_numOwnedNodes;
838  std::vector<LO> hi_to_lo_map;
839 
840  // Degree-1 variables
841  std::vector<size_t> lo_node_in_hi;
842 
843  // Degree-n variables
844  FCi lo_elemToHiRepresentativeNode;
845 
846  // Get Dirichlet unknown information
847  RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> > hi_isDirichletRow, hi_isDirichletCol;
848  Utilities::FindDirichletRowsAndPropagateToCols(A,hi_isDirichletRow, hi_isDirichletCol);
849 
850 #if 0
851  printf("[%d] isDirichletRow = ",A->getRowMap()->getComm()->getRank());
852  for(size_t i=0;i<hi_isDirichletRow->getMap()->getNodeNumElements(); i++)
853  printf("%d ",hi_isDirichletRow->getData(0)[i]);
854  printf("\n");
855  printf("[%d] isDirichletCol = ",A->getRowMap()->getComm()->getRank());
856  for(size_t i=0;i<hi_isDirichletCol->getMap()->getNodeNumElements(); i++)
857  printf("%d ",hi_isDirichletCol->getData(0)[i]);
858  printf("\n");
859  fflush(stdout);
860 #endif
861 
862  /*******************/
863  if(lo_degree == 1) {
864  // Get reference coordinates and the lo-to-hi injection list for the reference element
865  MueLuIntrepid::IntrepidGetP1NodeInHi(hi_basis,lo_node_in_hi,hi_DofCoords);
866 
867  // Generate lower-order element-to-node map
868  MueLuIntrepid::BuildLoElemToNode(*Pn_elemToNode,Pn_nodeIsOwned,lo_node_in_hi,hi_isDirichletCol->getData(0),*P1_elemToNode,P1_nodeIsOwned,hi_to_lo_map,P1_numOwnedNodes);
869  assert(hi_to_lo_map.size() == colMap->getNodeNumElements());
870  }
871  else {
872  // Get lo-order candidates
873  double threshold = 1e-10;
874  std::vector<std::vector<size_t> > candidates;
875  Kokkos::resize(hi_DofCoords,hi_basis->getCardinality(),hi_basis->getBaseCellTopology().getDimension());
876  hi_basis->getDofCoords(hi_DofCoords);
877 
878  MueLu::MueLuIntrepid::GenerateRepresentativeBasisNodes<Basis,FC>(*lo_basis,hi_DofCoords,threshold,candidates);
879 
880  // Generate the representative nodes
881  MueLu::MueLuIntrepid::GenerateLoNodeInHiViaGIDs(candidates,*Pn_elemToNode,colMap,lo_elemToHiRepresentativeNode);
882  MueLu::MueLuIntrepid::BuildLoElemToNodeViaRepresentatives(*Pn_elemToNode,Pn_nodeIsOwned,lo_elemToHiRepresentativeNode,*P1_elemToNode,P1_nodeIsOwned,hi_to_lo_map,P1_numOwnedNodes);
883  }
884  MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(coarseLevel,"pcoarsen: element to node map",P1_elemToNode);
885 
886  /*******************/
887  // Generate the P1_domainMap
888  // HOW: Since we know how many each proc has, we can use the non-uniform contiguous map constructor to do the work for us
889  RCP<const Map> P1_domainMap = MapFactory::Build(rowMap->lib(),Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),P1_numOwnedNodes,rowMap->getIndexBase(),rowMap->getComm());
890  MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(coarseLevel,"CoarseMap",P1_domainMap);
891 
892  // Generate the P1_columnMap
893  RCP<const Map> P1_colMap;
894  if(NumProc==1) P1_colMap = P1_domainMap;
895  else MueLuIntrepid::GenerateColMapFromImport<LO,GO,NO>(*Acrs.getCrsGraph()->getImporter(),hi_to_lo_map,*P1_domainMap,P1_nodeIsOwned.size(),P1_colMap);
896 
897  /*******************/
898  // Generate the coarsening
899  if(lo_degree == 1)
900  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);
901  else
902  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);
903 
904  /*******************/
905  // Zero out the Dirichlet rows in P
906  Utilities::ZeroDirichletRows(finalP,A_dirichletRows);
907 
908  /*******************/
909  // Build the nullspace
910  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P1_domainMap, fineNullspace->getNumVectors());
911  finalP->apply(*fineNullspace,*coarseNullspace,Teuchos::TRANS);
912  Set(coarseLevel, "Nullspace", coarseNullspace);
913 
914  // Level Set
915  if (!restrictionMode_) {
916  // The factory is in prolongation mode
917  Set(coarseLevel, "P", finalP);
918 
919  APparams->set("graph", finalP);
920  MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(coarseLevel,"AP reuse data",APparams);
921 
922  if (IsPrint(Statistics1)) {
923  RCP<ParameterList> params = rcp(new ParameterList());
924  params->set("printLoadBalancingInfo", true);
925  params->set("printCommInfo", true);
926  GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*finalP, "P", params);
927  }
928  } else {
929  // The factory is in restriction mode
930  RCP<Matrix> R;
931  {
932  SubFactoryMonitor m2(*this, "Transpose P", coarseLevel);
933  R = Utilities::Transpose(*finalP, true);
934  }
935 
936  Set(coarseLevel, "R", R);
937 
938  if (IsPrint(Statistics2)) {
939  RCP<ParameterList> params = rcp(new ParameterList());
940  params->set("printLoadBalancingInfo", true);
941  params->set("printCommInfo", true);
942  GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*R, "R", params);
943  }
944  }
945 
946  } //Build()
947 
948 } //namespace MueLu
949 
950 #endif // MUELU_IPCFACTORY_DEF_HPP
951 
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)
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
Teuchos::RCP< Intrepid2::Basis< KokkosExecutionSpace, Scalar, Scalar > > BasisFactory(const std::string &name, int &degree)
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)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
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
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
One-liner description of what is happening.
std::string tolower(const std::string &str)
Namespace for MueLu classes and methods.
MueLu::DefaultNode Node
Print even more statistics.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Kokkos::DynRankView< double, typename Node::device_type > SCFieldContainer
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)
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
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 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)
#define MUELU_LEVEL_SET_IF_REQUESTED_OR_KEPT(level, ename, entry)
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
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)
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
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
Kokkos::DynRankView< LocalOrdinal, typename Node::device_type > LOFieldContainer
Intrepid2::Basis< typename Node::device_type::execution_space, double, double > Basis
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 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 GenerateRepresentativeBasisNodes(const Basis &basis, const SCFieldContainer &ReferenceNodeLocations, const double threshold, std::vector< std::vector< size_t > > &representative_node_candidates)
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
#define SET_VALID_ENTRY(name)