MueLu  Version of the Day
MueLu_VariableDofLaplacianFactory_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 // Tobias Wiesner (tawiesn@sandia.gov)
42 // Ray Tuminaro (rstumin@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef PACKAGES_MUELU_SRC_GRAPH_MUELU_VARIABLEDOFLAPLACIANFACTORY_DEF_HPP_
48 #define PACKAGES_MUELU_SRC_GRAPH_MUELU_VARIABLEDOFLAPLACIANFACTORY_DEF_HPP_
49 
50 
51 #include "MueLu_Monitor.hpp"
52 
54 
55 namespace MueLu {
56 
57  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
59  RCP<ParameterList> validParamList = rcp(new ParameterList());
60 
61  validParamList->set< double > ("Advanced Dirichlet: threshold", 1e-5, "Drop tolerance for Dirichlet detection");
62  validParamList->set< double > ("Variable DOF amalgamation: threshold", 1.8e-9, "Drop tolerance for amalgamation process");
63  validParamList->set< int > ("maxDofPerNode", 1, "Maximum number of DOFs per node");
64 
65  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
66  validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for Coordinates");
67 
68  return validParamList;
69  }
70 
71  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
73 
74  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
76  Input(currentLevel, "A");
77  Input(currentLevel, "Coordinates");
78 
79  //if (currentLevel.GetLevelID() == 0) // TODO check for finest level (special treatment)
80  if (currentLevel.IsAvailable("DofPresent", NoFactory::get())) {
81  currentLevel.DeclareInput("DofPresent", NoFactory::get(), this);
82  }
83  }
84 
85  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
87  FactoryMonitor m(*this, "Build", currentLevel);
88  typedef Teuchos::ScalarTraits<SC> STS;
89 
90  const ParameterList & pL = GetParameterList();
91 
92  RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel, "A");
93 
94  Teuchos::RCP< const Teuchos::Comm< int > > comm = A->getRowMap()->getComm();
95  Xpetra::UnderlyingLib lib = A->getRowMap()->lib();
96 
97  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> dxMV;
98  RCP<dxMV> Coords = Get< RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > >(currentLevel, "Coordinates");
99 
100  int maxDofPerNode = pL.get<int>("maxDofPerNode");
101  Scalar dirDropTol = Teuchos::as<Scalar>(pL.get<double>("Advanced Dirichlet: threshold")); // "ML advnaced Dirichlet: threshold"
102  Scalar amalgDropTol = Teuchos::as<Scalar>(pL.get<double>("Variable DOF amalgamation: threshold")); //"variable DOF amalgamation: threshold")
103 
104  bool bHasZeroDiagonal = false;
105  Teuchos::ArrayRCP<const bool> dirOrNot = MueLu::Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DetectDirichletRowsExt(*A,bHasZeroDiagonal,STS::magnitude(dirDropTol));
106 
107  // check availability of DofPresent array
108  Teuchos::ArrayRCP<LocalOrdinal> dofPresent;
109  if (currentLevel.IsAvailable("DofPresent", NoFactory::get())) {
110  dofPresent = currentLevel.Get< Teuchos::ArrayRCP<LocalOrdinal> >("DofPresent", NoFactory::get());
111  } else {
112  // TAW: not sure about size of array. We cannot determine the expected size in the non-padded case correctly...
113  dofPresent = Teuchos::ArrayRCP<LocalOrdinal>(A->getRowMap()->getNodeNumElements(),1);
114  }
115 
116  // map[k] indicates that the kth dof in the variable dof matrix A would
117  // correspond to the map[k]th dof in the padded system. If, i.e., it is
118  // map[35] = 39 then dof no 35 in the variable dof matrix A corresponds to
119  // row map id 39 in an imaginary padded matrix Apadded.
120  // The padded system is never built but would be the associated matrix if
121  // every node had maxDofPerNode dofs.
122  std::vector<LocalOrdinal> map(A->getNodeNumRows());
123  this->buildPaddedMap(dofPresent, map, A->getNodeNumRows());
124 
125  // map of size of number of DOFs containing local node id (dof id -> node id, inclusive ghosted dofs/nodes)
126  std::vector<LocalOrdinal> myLocalNodeIds(A->getColMap()->getNodeNumElements()); // possible maximum (we need the ghost nodes, too)
127 
128  // assign the local node ids for the ghosted nodes
129  size_t nLocalNodes, nLocalPlusGhostNodes;
130  this->assignGhostLocalNodeIds(A->getRowMap(), A->getColMap(), myLocalNodeIds, map, maxDofPerNode, nLocalNodes, nLocalPlusGhostNodes, comm);
131 
132  //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)," ",0,false,10,false, true);
133 
134  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(dofPresent.size()) != Teuchos::as<size_t>(nLocalNodes * maxDofPerNode),MueLu::Exceptions::RuntimeError,"VariableDofLaplacianFactory: size of provided DofPresent array is " << dofPresent.size() << " but should be " << nLocalNodes * maxDofPerNode << " on the current processor.");
135 
136  // put content of assignGhostLocalNodeIds here...
137 
138  // fill nodal maps
139 
140  Teuchos::ArrayView< const GlobalOrdinal > myGids = A->getColMap()->getNodeElementList();
141 
142  // vector containing row/col gids of amalgamated matrix (with holes)
143 
144  size_t nLocalDofs = A->getRowMap()->getNodeNumElements();
145  size_t nLocalPlusGhostDofs = A->getColMap()->getNodeNumElements();
146 
147  // myLocalNodeIds (dof -> node)
148 
149  Teuchos::Array<GlobalOrdinal> amalgRowMapGIDs(nLocalNodes);
150  Teuchos::Array<GlobalOrdinal> amalgColMapGIDs(nLocalPlusGhostNodes);
151 
152  // initialize
153  size_t count = 0;
154  if (nLocalDofs > 0) {
155  amalgRowMapGIDs[count] = myGids[0];
156  amalgColMapGIDs[count] = myGids[0];
157  count++;
158  }
159 
160  for(size_t i = 1; i < nLocalDofs; i++) {
161  if (myLocalNodeIds[i] != myLocalNodeIds[i-1]) {
162  amalgRowMapGIDs[count] = myGids[i];
163  amalgColMapGIDs[count] = myGids[i];
164  count++;
165  }
166  }
167 
168  RCP<GOVector> tempAmalgColVec = GOVectorFactory::Build(A->getDomainMap());
169  {
170  Teuchos::ArrayRCP<GlobalOrdinal> tempAmalgColVecData = tempAmalgColVec->getDataNonConst(0);
171  for (size_t i = 0; i < A->getDomainMap()->getNodeNumElements(); i++)
172  tempAmalgColVecData[i] = amalgColMapGIDs[ myLocalNodeIds[i]];
173  }
174 
175  RCP<GOVector> tempAmalgColVecTarget = GOVectorFactory::Build(A->getColMap());
176  Teuchos::RCP<Import> dofImporter = ImportFactory::Build(A->getDomainMap(), A->getColMap());
177  tempAmalgColVecTarget->doImport(*tempAmalgColVec, *dofImporter, Xpetra::INSERT);
178 
179  {
180  Teuchos::ArrayRCP<const GlobalOrdinal> tempAmalgColVecBData = tempAmalgColVecTarget->getData(0);
181  // copy from dof vector to nodal vector
182  for (size_t i = 0; i < myLocalNodeIds.size(); i++)
183  amalgColMapGIDs[ myLocalNodeIds[i]] = tempAmalgColVecBData[i];
184  }
185 
186  Teuchos::RCP<Map> amalgRowMap = MapFactory::Build(lib,
187  Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),
188  amalgRowMapGIDs(), //View,
189  A->getRowMap()->getIndexBase(),
190  comm);
191 
192  Teuchos::RCP<Map> amalgColMap = MapFactory::Build(lib,
193  Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),
194  amalgColMapGIDs(), //View,
195  A->getRangeMap()->getIndexBase(),
196  comm);
197 
198  // end fill nodal maps
199 
200 
201  // start variable dof amalgamation
202 
203  Teuchos::RCP<CrsMatrixWrap> Awrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(A);
204  Teuchos::RCP<CrsMatrix> Acrs = Awrap->getCrsMatrix();
205  //Acrs->describe(*fancy, Teuchos::VERB_EXTREME);
206 
207  size_t nNonZeros = 0;
208  std::vector<bool> isNonZero(nLocalPlusGhostDofs,false);
209  std::vector<size_t> nonZeroList(nLocalPlusGhostDofs); // ???
210 
211  // also used in DetectDirichletExt
212  Teuchos::RCP<Vector> diagVecUnique = VectorFactory::Build(A->getRowMap());
213  Teuchos::RCP<Vector> diagVec = VectorFactory::Build(A->getColMap());
214  A->getLocalDiagCopy(*diagVecUnique);
215  diagVec->doImport(*diagVecUnique, *dofImporter, Xpetra::INSERT);
216  Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
217 
218  Teuchos::ArrayRCP<const size_t> rowptr(Acrs->getNodeNumRows());
219  Teuchos::ArrayRCP<const LocalOrdinal> colind(Acrs->getNodeNumEntries());
220  Teuchos::ArrayRCP<const Scalar> values(Acrs->getNodeNumEntries());
221  Acrs->getAllValues(rowptr, colind, values);
222 
223 
224  // create arrays for amalgamated matrix
225  Teuchos::ArrayRCP<size_t> amalgRowPtr(nLocalNodes+1);
226  Teuchos::ArrayRCP<LocalOrdinal> amalgCols(rowptr[rowptr.size()-1]);
227 
228  LocalOrdinal oldBlockRow = 0;
229  LocalOrdinal blockRow = 0;
230  LocalOrdinal blockColumn = 0;
231 
232  size_t newNzs = 0;
233  amalgRowPtr[0] = newNzs;
234 
235  bool doNotDrop = false;
236  if (amalgDropTol == Teuchos::ScalarTraits<Scalar>::zero()) doNotDrop = true;
237  if (values.size() == 0) doNotDrop = true;
238 
239  for(decltype(rowptr.size()) i = 0; i < rowptr.size()-1; i++) {
240  blockRow = std::floor<LocalOrdinal>( map[i] / maxDofPerNode);
241  if (blockRow != oldBlockRow) {
242  // zero out info recording nonzeros in oldBlockRow
243  for(size_t j = 0; j < nNonZeros; j++) isNonZero[nonZeroList[j]] = false;
244  nNonZeros = 0;
245  amalgRowPtr[blockRow] = newNzs; // record start of next row
246  }
247  for (size_t j = rowptr[i]; j < rowptr[i+1]; j++) {
248  if(doNotDrop == true ||
249  ( STS::magnitude(values[j] / STS::magnitude(sqrt(STS::magnitude(diagVecData[i]) * STS::magnitude(diagVecData[colind[j]]))) ) >= STS::magnitude(amalgDropTol) )) {
250  blockColumn = myLocalNodeIds[colind[j]];
251  if(isNonZero[blockColumn] == false) {
252  isNonZero[blockColumn] = true;
253  nonZeroList[nNonZeros++] = blockColumn;
254  amalgCols[newNzs++] = blockColumn;
255  }
256  }
257  }
258  oldBlockRow = blockRow;
259  }
260  amalgRowPtr[blockRow+1] = newNzs;
261 
262  TEUCHOS_TEST_FOR_EXCEPTION((blockRow+1 != Teuchos::as<LO>(nLocalNodes)) && (nLocalNodes !=0), MueLu::Exceptions::RuntimeError, "VariableDofsPerNodeAmalgamation: error, computed # block rows (" << blockRow+1 <<") != nLocalNodes (" << nLocalNodes <<")");
263 
264  amalgCols.resize(amalgRowPtr[nLocalNodes]);
265 
266  // end variableDofAmalg
267 
268  // begin rm differentDofsCrossings
269 
270  // Remove matrix entries (i,j) where the ith node and the jth node have
271  // different dofs that are 'present'
272  // Specifically, on input:
273  // dofPresent[i*maxDofPerNode+k] indicates whether or not the kth
274  // dof at the ith node is present in the
275  // variable dof matrix (e.g., the ith node
276  // has an air pressure dof). true means
277  // the dof is present while false means it
278  // is not.
279  // We create a unique id for the ith node (i.e. uniqueId[i]) via
280  // sum_{k=0 to maxDofPerNode-1} dofPresent[i*maxDofPerNode+k]*2^k
281  // and use this unique idea to remove entries (i,j) when uniqueId[i]!=uniqueId[j]
282 
283  Teuchos::ArrayRCP<LocalOrdinal> uniqueId(nLocalPlusGhostNodes); // unique id associated with DOF
284  std::vector<bool> keep(amalgRowPtr[amalgRowPtr.size()-1],true); // keep connection associated with node
285 
286  size_t ii = 0; // iteration index for present dofs
287  for(decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size()-1; i++) {
288  LocalOrdinal temp = 1; // basis for dof-id
289  uniqueId[i] = 0;
290  for (decltype(maxDofPerNode) j = 0; j < maxDofPerNode; j++) {
291  if (dofPresent[ii++]) uniqueId[i] += temp; // encode dof to be present
292  temp = temp * 2; // check next dof
293  }
294  }
295 
296  Teuchos::RCP<Import> nodeImporter = ImportFactory::Build(amalgRowMap, amalgColMap);
297 
298  RCP<LOVector> nodeIdSrc = Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(amalgRowMap,true);
299  RCP<LOVector> nodeIdTarget = Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(amalgColMap,true);
300 
301  Teuchos::ArrayRCP< LocalOrdinal > nodeIdSrcData = nodeIdSrc->getDataNonConst(0);
302  for(decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size()-1; i++) {
303  nodeIdSrcData[i] = uniqueId[i];
304  }
305 
306  nodeIdTarget->doImport(*nodeIdSrc, *nodeImporter, Xpetra::INSERT);
307 
308  Teuchos::ArrayRCP< const LocalOrdinal > nodeIdTargetData = nodeIdTarget->getData(0);
309  for(decltype(uniqueId.size()) i = 0; i < uniqueId.size(); i++) {
310  uniqueId[i] = nodeIdTargetData[i];
311  }
312 
313  // nodal comm uniqueId, myLocalNodeIds
314 
315  // uniqueId now should contain ghosted data
316 
317  for(decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size()-1; i++) {
318  for(size_t j = amalgRowPtr[i]; j < amalgRowPtr[i+1]; j++) {
319  if (uniqueId[i] != uniqueId[amalgCols[j]]) keep [j] = false;
320  }
321  }
322 
323  // squeeze out hard-coded zeros from CSR arrays
324  Teuchos::ArrayRCP<Scalar> amalgVals;
325  this->squeezeOutNnzs(amalgRowPtr,amalgCols,amalgVals,keep);
326 
327  typedef Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> dxMVf;
328  RCP<dxMV> ghostedCoords = dxMVf::Build(amalgColMap,Coords->getNumVectors());
329 
330  TEUCHOS_TEST_FOR_EXCEPTION(amalgRowMap->getNodeNumElements() != Coords->getMap()->getNodeNumElements(), MueLu::Exceptions::RuntimeError, "MueLu::VariableDofLaplacianFactory: the number of Coordinates and amalgamated nodes is inconsistent.");
331 
332  // Coords might live on a special nodeMap with consecutive ids (the natural numbering)
333  // The amalgRowMap might have the same number of entries, but with holes in the ids.
334  // e.g. 0,3,6,9,... as GIDs.
335  // We need the ghosted Coordinates in the buildLaplacian routine. But we access the data
336  // through getData only, i.e., the global ids are not interesting as long as we do not change
337  // the ordering of the entries
338  Coords->replaceMap(amalgRowMap);
339  ghostedCoords->doImport(*Coords, *nodeImporter, Xpetra::INSERT);
340 
341  Teuchos::ArrayRCP<Scalar> lapVals(amalgRowPtr[nLocalNodes]);
342  this->buildLaplacian(amalgRowPtr, amalgCols, lapVals, Coords->getNumVectors(), ghostedCoords);
343 
344  // sort column GIDs
345  for(decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size()-1; i++) {
346  size_t j = amalgRowPtr[i];
347  this->MueLu_az_sort<LocalOrdinal>(&(amalgCols[j]), amalgRowPtr[i+1] - j, NULL, &(lapVals[j]));
348  }
349 
350  // Caluclate status array for next level
351  Teuchos::Array<char> status(nLocalNodes * maxDofPerNode);
352 
353  // dir or not Teuchos::ArrayRCP<const bool> dirOrNot
354  for(decltype(status.size()) i = 0; i < status.size(); i++) status[i] = 's';
355  for(decltype(status.size()) i = 0; i < status.size(); i++) {
356  if(dofPresent[i] == false) status[i] = 'p';
357  }
358  if(dirOrNot.size() > 0) {
359  for(decltype(map.size()) i = 0; i < map.size(); i++) {
360  if(dirOrNot[i] == true){
361  status[map[i]] = 'd';
362  }
363  }
364  }
365  Set(currentLevel,"DofStatus",status);
366 
367  // end status array
368 
369  Teuchos::RCP<CrsMatrix> lapCrsMat = CrsMatrixFactory::Build(amalgRowMap, amalgColMap, 10); // TODO better approx for max nnz per row
370 
371  for (size_t i = 0; i < nLocalNodes; i++) {
372  lapCrsMat->insertLocalValues(i, amalgCols.view(amalgRowPtr[i],amalgRowPtr[i+1]-amalgRowPtr[i]),
373  lapVals.view(amalgRowPtr[i],amalgRowPtr[i+1]-amalgRowPtr[i]));
374  }
375  lapCrsMat->fillComplete(amalgRowMap,amalgRowMap);
376 
377  //lapCrsMat->describe(*fancy, Teuchos::VERB_EXTREME);
378 
379  Teuchos::RCP<Matrix> lapMat = Teuchos::rcp(new CrsMatrixWrap(lapCrsMat));
380  Set(currentLevel,"A",lapMat);
381  }
382 
383  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
384  void VariableDofLaplacianFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::buildLaplacian(const Teuchos::ArrayRCP<size_t>& rowPtr, const Teuchos::ArrayRCP<LocalOrdinal>& cols, Teuchos::ArrayRCP<Scalar>& vals,const size_t& numdim, const RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LocalOrdinal,GlobalOrdinal,Node> > & ghostedCoords) const {
385  TEUCHOS_TEST_FOR_EXCEPTION(numdim != 2 && numdim !=3, MueLu::Exceptions::RuntimeError,"buildLaplacian only works for 2d or 3d examples. numdim = " << numdim);
386 
387  if(numdim == 2) { // 2d
388  Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits<Scalar>::magnitudeType > x = ghostedCoords->getData(0);
389  Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits<Scalar>::magnitudeType > y = ghostedCoords->getData(1);
390 
391  for(decltype(rowPtr.size()) i = 0; i < rowPtr.size() - 1; i++) {
392  Scalar sum = Teuchos::ScalarTraits<Scalar>::zero();
393  LocalOrdinal diag = -1;
394  for(size_t j = rowPtr[i]; j < rowPtr[i+1]; j++) {
395  if(cols[j] != Teuchos::as<LO>(i)){
396  vals[j] = std::sqrt( (x[i]-x[cols[j]]) * (x[i]-x[cols[j]]) +
397  (y[i]-y[cols[j]]) * (y[i]-y[cols[j]]) );
398  TEUCHOS_TEST_FOR_EXCEPTION(vals[j] == Teuchos::ScalarTraits<Scalar>::zero(), MueLu::Exceptions::RuntimeError, "buildLaplacian: error, " << i << " and " << cols[j] << " have same coordinates: " << x[i] << " and " << y[i]);
399  vals[j] = -Teuchos::ScalarTraits<SC>::one()/vals[j];
400  sum = sum - vals[j];
401  }
402  else diag = j;
403  }
404  if(sum == Teuchos::ScalarTraits<Scalar>::zero()) sum = Teuchos::ScalarTraits<Scalar>::one();
405  TEUCHOS_TEST_FOR_EXCEPTION(diag == -1, MueLu::Exceptions::RuntimeError, "buildLaplacian: error, row " << i << " has zero diagonal!");
406 
407  vals[diag] = sum;
408  }
409  } else { // 3d
410  Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits<Scalar>::magnitudeType > x = ghostedCoords->getData(0);
411  Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits<Scalar>::magnitudeType > y = ghostedCoords->getData(1);
412  Teuchos::ArrayRCP< const typename Teuchos::ScalarTraits<Scalar>::magnitudeType > z = ghostedCoords->getData(2);
413 
414  for(decltype(rowPtr.size()) i = 0; i < rowPtr.size() - 1; i++) {
415  Scalar sum = Teuchos::ScalarTraits<Scalar>::zero();
416  LocalOrdinal diag = -1;
417  for(size_t j = rowPtr[i]; j < rowPtr[i+1]; j++) {
418  if(cols[j] != Teuchos::as<LO>(i)){
419  vals[j] = std::sqrt( (x[i]-x[cols[j]]) * (x[i]-x[cols[j]]) +
420  (y[i]-y[cols[j]]) * (y[i]-y[cols[j]]) +
421  (z[i]-z[cols[j]]) * (z[i]-z[cols[j]]) );
422 
423  TEUCHOS_TEST_FOR_EXCEPTION(vals[j] == Teuchos::ScalarTraits<Scalar>::zero(), MueLu::Exceptions::RuntimeError, "buildLaplacian: error, " << i << " and " << cols[j] << " have same coordinates: " << x[i] << " and " << y[i] << " and " << z[i]);
424 
425  vals[j] = -Teuchos::ScalarTraits<SC>::one()/vals[j];
426  sum = sum - vals[j];
427  }
428  else diag = j;
429  }
430  if(sum == Teuchos::ScalarTraits<Scalar>::zero()) sum = Teuchos::ScalarTraits<Scalar>::one();
431  TEUCHOS_TEST_FOR_EXCEPTION(diag == -1, MueLu::Exceptions::RuntimeError, "buildLaplacian: error, row " << i << " has zero diagonal!");
432 
433  vals[diag] = sum;
434  }
435  }
436  }
437 
438  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
439  void VariableDofLaplacianFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::squeezeOutNnzs(Teuchos::ArrayRCP<size_t> & rowPtr, Teuchos::ArrayRCP<LocalOrdinal> & cols, Teuchos::ArrayRCP<Scalar> & vals, const std::vector<bool>& keep) const {
440  // get rid of nonzero entries that have 0's in them and properly change
441  // the row ptr array to reflect this removal (either vals == NULL or vals != NULL)
442  // Note, the arrays are squeezed. No memory is freed.
443 
444  size_t count = 0;
445 
446  size_t nRows = rowPtr.size()-1;
447  if(vals.size() > 0) {
448  for(size_t i = 0; i < nRows; i++) {
449  size_t newStart = count;
450  for(size_t j = rowPtr[i]; j < rowPtr[i+1]; j++) {
451  if(vals[j] != Teuchos::ScalarTraits<Scalar>::zero()) {
452  cols[count ] = cols[j];
453  vals[count++] = vals[j];
454  }
455  }
456  rowPtr[i] = newStart;
457  }
458  } else {
459  for (size_t i = 0; i < nRows; i++) {
460  size_t newStart = count;
461  for(size_t j = rowPtr[i]; j < rowPtr[i+1]; j++) {
462  if (keep[j] == true) {
463  cols[count++] = cols[j];
464  }
465  }
466  rowPtr[i] = newStart;
467  }
468  }
469  rowPtr[nRows] = count;
470  }
471 
472  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
473  void VariableDofLaplacianFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::buildPaddedMap(const Teuchos::ArrayRCP<const LocalOrdinal> & dofPresent, std::vector<LocalOrdinal> & map, size_t nDofs) const {
474  size_t count = 0;
475  for (decltype(dofPresent.size()) i = 0; i < dofPresent.size(); i++)
476  if(dofPresent[i] == 1) map[count++] = Teuchos::as<LocalOrdinal>(i);
477  TEUCHOS_TEST_FOR_EXCEPTION(nDofs != count, MueLu::Exceptions::RuntimeError, "VariableDofLaplacianFactory::buildPaddedMap: #dofs in dofPresent does not match the expected value (number of rows of A): " << nDofs << " vs. " << count);
478  }
479 
480  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
481  void VariableDofLaplacianFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::assignGhostLocalNodeIds(const Teuchos::RCP<const Map> & rowDofMap, const Teuchos::RCP<const Map> & colDofMap, std::vector<LocalOrdinal> & myLocalNodeIds, const std::vector<LocalOrdinal> & dofMap, size_t maxDofPerNode, size_t& nLocalNodes, size_t& nLocalPlusGhostNodes, Teuchos::RCP< const Teuchos::Comm< int > > comm) const {
482 
483  size_t nLocalDofs = rowDofMap->getNodeNumElements();
484  size_t nLocalPlusGhostDofs = colDofMap->getNodeNumElements(); // TODO remove parameters
485 
486  // create importer for dof-based information
487  Teuchos::RCP<Import> importer = ImportFactory::Build(rowDofMap, colDofMap);
488 
489  // create a vector living on column map of A (dof based)
490  Teuchos::RCP<LOVector> localNodeIdsTemp = LOVectorFactory::Build(rowDofMap,true);
491  Teuchos::RCP<LOVector> localNodeIds = LOVectorFactory::Build(colDofMap,true);
492 
493  // fill local dofs (padded local ids)
494  {
495  Teuchos::ArrayRCP< LocalOrdinal > localNodeIdsTempData = localNodeIdsTemp->getDataNonConst(0);
496  for(size_t i = 0; i < localNodeIdsTemp->getLocalLength(); i++)
497  localNodeIdsTempData[i] = std::floor<LocalOrdinal>( dofMap[i] / maxDofPerNode );
498  }
499 
500  localNodeIds->doImport(*localNodeIdsTemp, *importer, Xpetra::INSERT);
501  Teuchos::ArrayRCP< const LocalOrdinal > localNodeIdsData = localNodeIds->getData(0);
502 
503  // Note: localNodeIds contains local ids for the padded version as vector values
504 
505 
506  // we use Scalar instead of int as type
507  Teuchos::RCP<LOVector> myProcTemp = LOVectorFactory::Build(rowDofMap,true);
508  Teuchos::RCP<LOVector> myProc = LOVectorFactory::Build(colDofMap,true);
509 
510  // fill local dofs (padded local ids)
511  {
512  Teuchos::ArrayRCP< LocalOrdinal > myProcTempData = myProcTemp->getDataNonConst(0);
513  for(size_t i = 0; i < myProcTemp->getLocalLength(); i++)
514  myProcTempData[i] = Teuchos::as<LocalOrdinal>(comm->getRank());
515  }
516  myProc->doImport(*myProcTemp, *importer, Xpetra::INSERT);
517  Teuchos::ArrayRCP<LocalOrdinal> myProcData = myProc->getDataNonConst(0); // we have to modify the data (therefore the non-const version)
518 
519  // At this point, the ghost part of localNodeIds corresponds to the local ids
520  // associated with the current owning processor. We want to convert these to
521  // local ids associated with the processor on which these are ghosts.
522  // Thus we have to re-number them. In doing this re-numbering we must make sure
523  // that we find all ghosts with the same id & proc and assign a unique local
524  // id to this group (id&proc). To do this find, we sort all ghost entries in
525  // localNodeIds that are owned by the same processor. Then we can look for
526  // duplicates (i.e., several ghost entries corresponding to dofs with the same
527  // node id) easily and make sure these are all assigned to the same local id.
528  // To do the sorting we'll make a temporary copy of the ghosts via tempId and
529  // tempProc and sort this multiple times for each group owned by the same proc.
530 
531 
532  std::vector<size_t> location(nLocalPlusGhostDofs - nLocalDofs + 1);
533  std::vector<size_t> tempId (nLocalPlusGhostDofs - nLocalDofs + 1);
534  std::vector<size_t> tempProc(nLocalPlusGhostDofs - nLocalDofs + 1);
535 
536  size_t notProcessed = nLocalDofs; // iteration index over all ghosted dofs
537  size_t tempIndex = 0;
538  size_t first = tempIndex;
539  LocalOrdinal neighbor;
540 
541  while (notProcessed < nLocalPlusGhostDofs) {
542  neighbor = myProcData[notProcessed]; // get processor id of not-processed element
543  first = tempIndex;
544  location[tempIndex] = notProcessed;
545  tempId[tempIndex++] = localNodeIdsData[notProcessed];
546  myProcData[notProcessed] = -1 - neighbor;
547 
548  for(size_t i = notProcessed + 1; i < nLocalPlusGhostDofs; i++) {
549  if(myProcData[i] == neighbor) {
550  location[tempIndex] = i;
551  tempId[tempIndex++] = localNodeIdsData[i];
552  myProcData[i] = -1; // mark as visited
553  }
554  }
555  this->MueLu_az_sort<size_t>(&(tempId[first]), tempIndex - first, &(location[first]), NULL);
556  for(size_t i = first; i < tempIndex; i++) tempProc[i] = neighbor;
557 
558  // increment index. Find next notProcessed dof index corresponding to first non-visited element
559  notProcessed++;
560  while ( (notProcessed < nLocalPlusGhostDofs) && (myProcData[notProcessed] < 0))
561  notProcessed++;
562  }
563  TEUCHOS_TEST_FOR_EXCEPTION(tempIndex != nLocalPlusGhostDofs-nLocalDofs, MueLu::Exceptions::RuntimeError,"Number of nonzero ghosts is inconsistent.");
564 
565  // Now assign ids to all ghost nodes (giving the same id to those with the
566  // same myProc[] and the same local id on the proc that actually owns the
567  // variable associated with the ghost
568 
569  nLocalNodes = 0; // initialize return value
570  if(nLocalDofs > 0) nLocalNodes = localNodeIdsData[nLocalDofs-1] + 1;
571 
572  nLocalPlusGhostNodes = nLocalNodes; // initialize return value
573  if(nLocalDofs < nLocalPlusGhostDofs) nLocalPlusGhostNodes++; // 1st ghost node is unique (not accounted for). number will be increased later, if there are more ghost nodes
574 
575  // check if two adjacent ghost dofs correspond to different nodes. To do this,
576  // check if they are from different processors or whether they have different
577  // local node ids
578 
579  // loop over all (remaining) ghost dofs
580  for (size_t i = nLocalDofs+1; i < nLocalPlusGhostDofs; i++) {
581  size_t lagged = nLocalPlusGhostNodes-1;
582 
583  // i is a new unique ghost node (not already accounted for)
584  if ((tempId[i-nLocalDofs] != tempId[i-1-nLocalDofs]) ||
585  (tempProc[i-nLocalDofs] != tempProc[i-1-nLocalDofs]))
586  nLocalPlusGhostNodes++; // update number of ghost nodes
587  tempId[i-1-nLocalDofs] = lagged;
588  }
589  if (nLocalPlusGhostDofs > nLocalDofs)
590  tempId[nLocalPlusGhostDofs-1-nLocalDofs] = nLocalPlusGhostNodes - 1;
591 
592  // fill myLocalNodeIds array. Start with local part (not ghosted)
593  for(size_t i = 0; i < nLocalDofs; i++)
594  myLocalNodeIds[i] = std::floor<LocalOrdinal>( dofMap[i] / maxDofPerNode );
595 
596  // copy ghosted nodal ids into myLocalNodeIds
597  for(size_t i = nLocalDofs; i < nLocalPlusGhostDofs; i++)
598  myLocalNodeIds[location[i-nLocalDofs]] = tempId[i-nLocalDofs];
599 
600  }
601 
602 } /* MueLu */
603 
604 
605 #endif /* PACKAGES_MUELU_SRC_GRAPH_MUELU_VARIABLEDOFLAPLACIANFACTORY_DEF_HPP_ */
MueLu::DefaultLocalOrdinal LocalOrdinal
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
void buildLaplacian(const Teuchos::ArrayRCP< size_t > &rowPtr, const Teuchos::ArrayRCP< LocalOrdinal > &cols, Teuchos::ArrayRCP< Scalar > &vals, const size_t &numdim, const RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > &ghostedCoords) const
MueLu::DefaultNode Node
static const NoFactory * get()
void assignGhostLocalNodeIds(const Teuchos::RCP< const Map > &rowDofMap, const Teuchos::RCP< const Map > &colDofMap, std::vector< LocalOrdinal > &myLocalNodeIds, const std::vector< LocalOrdinal > &dofMap, size_t maxDofPerNode, size_t &nLocalNodes, size_t &nLocalPlusGhostNodes, Teuchos::RCP< const Teuchos::Comm< int > > comm) const
void buildPaddedMap(const Teuchos::ArrayRCP< const LocalOrdinal > &dofPresent, std::vector< LocalOrdinal > &map, size_t nDofs) const
void squeezeOutNnzs(Teuchos::ArrayRCP< size_t > &rowPtr, Teuchos::ArrayRCP< LocalOrdinal > &cols, Teuchos::ArrayRCP< Scalar > &vals, const std::vector< bool > &keep) const
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void Build(Level &currentLevel) const
Build an object with this factory.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Exception throws to report errors in the internal logical of the program.
void DeclareInput(Level &currentLevel) const
Input.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()