MueLu Version of the Day
Loading...
Searching...
No Matches
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
55namespace 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()->getLocalNumElements(),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->getLocalNumRows());
123 this->buildPaddedMap(dofPresent, map, A->getLocalNumRows());
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()->getLocalNumElements()); // 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()->getLocalElementList();
141
142 // vector containing row/col gids of amalgamated matrix (with holes)
143
144 size_t nLocalDofs = A->getRowMap()->getLocalNumElements();
145 size_t nLocalPlusGhostDofs = A->getColMap()->getLocalNumElements();
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()->getLocalNumElements(); 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->getLocalNumRows());
219 Teuchos::ArrayRCP<const LocalOrdinal> colind(Acrs->getLocalNumEntries());
220 Teuchos::ArrayRCP<const Scalar> values(Acrs->getLocalNumEntries());
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->getLocalNumElements() != Coords->getMap()->getLocalNumElements(), 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->getLocalNumElements();
484 size_t nLocalPlusGhostDofs = colDofMap->getLocalNumElements(); // 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
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static const NoFactory * get()
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
void Build(Level &currentLevel) const
Build an object with this factory.
void squeezeOutNnzs(Teuchos::ArrayRCP< size_t > &rowPtr, Teuchos::ArrayRCP< LocalOrdinal > &cols, Teuchos::ArrayRCP< Scalar > &vals, const std::vector< bool > &keep) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &currentLevel) const
Input.
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
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
Namespace for MueLu classes and methods.