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