MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_LineDetectionFactory_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_LINEDETECTIONFACTORY_DEF_HPP
47#define MUELU_LINEDETECTIONFACTORY_DEF_HPP
48
49#include <Xpetra_Matrix.hpp>
50//#include <Xpetra_MatrixFactory.hpp>
51
53
54#include "MueLu_FactoryManager.hpp"
55#include "MueLu_Level.hpp"
56#include "MueLu_MasterList.hpp"
57#include "MueLu_Monitor.hpp"
58
59namespace MueLu {
60
61 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63 RCP<ParameterList> validParamList = rcp(new ParameterList());
64
65#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
66 SET_VALID_ENTRY("linedetection: orientation");
67 SET_VALID_ENTRY("linedetection: num layers");
68#undef SET_VALID_ENTRY
69
70 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
71 validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for coorindates");
72
73 return validParamList;
74 }
75
76 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
78 Input(currentLevel, "A");
79
80 // The factory needs the information about the number of z-layers. While this information is
81 // provided by the user for the finest level, the factory itself is responsible to provide the
82 // corresponding information on the coarser levels. Since a factory cannot be dependent on itself
83 // we use the NoFactory class as generator class, but remove the UserData keep flag, such that
84 // "NumZLayers" is part of the request/release mechanism.
85 // Please note, that this prevents us from having several (independent) CoarsePFactory instances!
86 // TODO: allow factory to dependent on self-generated data for TwoLevelFactories -> introduce ExpertRequest/Release in Level
87 currentLevel.DeclareInput("NumZLayers", NoFactory::get(), this);
88 currentLevel.RemoveKeepFlag("NumZLayers", NoFactory::get(), MueLu::UserData);
89 }
90
91 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93 FactoryMonitor m(*this, "Line detection (Ray style)", currentLevel);
94
95 LO NumZDir = 0;
96 RCP<CoordinateMultiVector> fineCoords;
97 ArrayRCP<coordinate_type> x, y, z;
98 coordinate_type *xptr = NULL, *yptr = NULL, *zptr = NULL;
99
100 // obtain general variables
101 RCP<Matrix> A = Get< RCP<Matrix> > (currentLevel, "A");
102 LO BlkSize = A->GetFixedBlockSize();
103 RCP<const Map> rowMap = A->getRowMap();
104 LO Ndofs = rowMap->getLocalNumElements();
105 LO Nnodes = Ndofs/BlkSize;
106
107 // collect information provided by user
108 const ParameterList& pL = GetParameterList();
109 const std::string lineOrientation = pL.get<std::string>("linedetection: orientation");
110
111 // interpret "line orientation" parameter provided by the user on the finest level
112 if(currentLevel.GetLevelID() == 0) {
113 if(lineOrientation=="vertical")
114 Zorientation_ = VERTICAL;
115 else if (lineOrientation=="horizontal")
116 Zorientation_ = HORIZONTAL;
117 else if (lineOrientation=="coordinates")
118 Zorientation_ = GRID_SUPPLIED;
119 else
120 TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "LineDetectionFactory: The parameter 'semicoarsen: line orientation' must be either 'vertical', 'horizontal' or 'coordinates'.");
121 }
122
123 //TEUCHOS_TEST_FOR_EXCEPTION(Zorientation_!=VERTICAL, Exceptions::RuntimeError, "LineDetectionFactory: The 'horizontal' or 'coordinates' have not been tested!!!. Please remove this exception check and carefully test these modes!");
124
125 // obtain number of z layers (variable over levels)
126 // This information is user-provided on the finest level and transferred to the coarser
127 // levels by the SemiCoarsenPFactor using the internal "NumZLayers" variable.
128 if(currentLevel.GetLevelID() == 0) {
129 if(currentLevel.IsAvailable("NumZLayers", NoFactory::get())) {
130 NumZDir = currentLevel.Get<LO>("NumZLayers", NoFactory::get()); //obtain info
131 GetOStream(Runtime1) << "Number of layers for line detection: " << NumZDir << " (information from Level(0))" << std::endl;
132 } else {
133 // check whether user provides information or it can be reconstructed from coordinates
134 NumZDir = pL.get<LO>("linedetection: num layers");
135 if(NumZDir == -1) {
136 bool CoordsAvail = currentLevel.IsAvailable("Coordinates");
137
138 if (CoordsAvail == true) {
139 // try to reconstruct the number of layers from coordinates
140 fineCoords = Get< RCP<CoordinateMultiVector> > (currentLevel, "Coordinates");
141 TEUCHOS_TEST_FOR_EXCEPTION(fineCoords->getNumVectors() != 3, Exceptions::RuntimeError, "Three coordinates arrays must be supplied if line detection orientation not given.");
142 x = fineCoords->getDataNonConst(0);
143 y = fineCoords->getDataNonConst(1);
144 z = fineCoords->getDataNonConst(2);
145 xptr = x.getRawPtr();
146 yptr = y.getRawPtr();
147 zptr = z.getRawPtr();
148
149 LO NumCoords = Ndofs/BlkSize;
150
151 /* sort coordinates so that we can order things according to lines */
152 Teuchos::ArrayRCP<LO> TOrigLoc= Teuchos::arcp<LO>(NumCoords); LO* OrigLoc= TOrigLoc.getRawPtr();
153 Teuchos::ArrayRCP<coordinate_type> Txtemp = Teuchos::arcp<coordinate_type>(NumCoords); coordinate_type* xtemp = Txtemp.getRawPtr();
154 Teuchos::ArrayRCP<coordinate_type> Tytemp = Teuchos::arcp<coordinate_type>(NumCoords); coordinate_type* ytemp = Tytemp.getRawPtr();
155 Teuchos::ArrayRCP<coordinate_type> Tztemp = Teuchos::arcp<coordinate_type>(NumCoords); coordinate_type* ztemp = Tztemp.getRawPtr();
156
157 // sort coordinates in {x,y,z}vals (returned in {x,y,z}temp) so that we can order things according to lines
158 // switch x and y coordinates for semi-coarsening...
159 sort_coordinates(NumCoords, OrigLoc, xptr, yptr, zptr, xtemp, ytemp, ztemp, true);
160
161 /* go through each vertical line and populate blockIndices so all */
162 /* dofs within a PDE within a vertical line correspond to one block.*/
163 LO NumBlocks = 0;
164 LO NumNodesPerVertLine = 0;
165 LO index = 0;
166
167 while ( index < NumCoords ) {
168 coordinate_type xfirst = xtemp[index]; coordinate_type yfirst = ytemp[index];
169 LO next = index+1;
170 while ( (next != NumCoords) && (xtemp[next] == xfirst) &&
171 (ytemp[next] == yfirst))
172 next++;
173 if (NumBlocks == 0) {
174 NumNodesPerVertLine = next-index;
175 }
176 // the number of vertical lines must be the same on all processors
177 // TAW: Sep 14 2015: or zero as we allow "empty" processors
178 //TEUCHOS_TEST_FOR_EXCEPTION(next-index != NumNodesPerVertLine,Exceptions::RuntimeError, "Error code only works for constant block size now!!!\n");
179 NumBlocks++;
180 index = next;
181 }
182
183 NumZDir = NumNodesPerVertLine;
184 GetOStream(Runtime1) << "Number of layers for line detection: " << NumZDir << " (information reconstructed from provided node coordinates)" << std::endl;
185 } else {
186 TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "LineDetectionFactory: BuildP: User has to provide valid number of layers (e.g. using the 'line detection: num layers' parameter).");
187 }
188 } else {
189 GetOStream(Runtime1) << "Number of layers for line detection: " << NumZDir << " (information provided by user through 'line detection: num layers')" << std::endl;
190 }
191 } // end else (user provides information or can be reconstructed) on finest level
192 } else {
193 // coarse level information
194 // TODO get rid of NoFactory here and use SemiCoarsenPFactory as source of NumZLayers instead.
195 if(currentLevel.IsAvailable("NumZLayers", NoFactory::get())) {
196 NumZDir = currentLevel.Get<LO>("NumZLayers", NoFactory::get()); //obtain info
197 GetOStream(Runtime1) << "Number of layers for line detection: " << NumZDir << std::endl;
198 } else {
199 TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "LineDetectionFactory: BuildP: No NumZLayers variable found. This cannot be.");
200 }
201 }
202
203 // plausibility check and further variable collection
204 if (Zorientation_ == GRID_SUPPLIED) { // On finest level, fetch user-provided coordinates if available...
205 bool CoordsAvail = currentLevel.IsAvailable("Coordinates");
206
207 if (CoordsAvail == false) {
208 if (currentLevel.GetLevelID() == 0)
209 throw Exceptions::RuntimeError("Coordinates must be supplied if line detection orientation not given.");
210 else
211 throw Exceptions::RuntimeError("Coordinates not generated by previous invocation of LineDetectionFactory's BuildP() method.");
212 }
213 fineCoords = Get< RCP<CoordinateMultiVector> > (currentLevel, "Coordinates");
214 TEUCHOS_TEST_FOR_EXCEPTION(fineCoords->getNumVectors() != 3, Exceptions::RuntimeError, "Three coordinates arrays must be supplied if line detection orientation not given.");
215 x = fineCoords->getDataNonConst(0);
216 y = fineCoords->getDataNonConst(1);
217 z = fineCoords->getDataNonConst(2);
218 xptr = x.getRawPtr();
219 yptr = y.getRawPtr();
220 zptr = z.getRawPtr();
221 }
222
223 // perform line detection
224 if (NumZDir > 0) {
225 LO *LayerId, *VertLineId;
226 Teuchos::ArrayRCP<LO> TLayerId = Teuchos::arcp<LO>(Nnodes); LayerId = TLayerId.getRawPtr();
227 Teuchos::ArrayRCP<LO> TVertLineId= Teuchos::arcp<LO>(Nnodes); VertLineId = TVertLineId.getRawPtr();
228
229 NumZDir = ML_compute_line_info(LayerId, VertLineId, Ndofs, BlkSize,
230 Zorientation_, NumZDir,xptr,yptr,zptr, *(rowMap->getComm()));
231 //it is NumZDir=NCLayers*NVertLines*DofsPerNode;
232
233 // store output data on current level
234 // The line detection data is used by the SemiCoarsenPFactory and the line smoothers in Ifpack/Ifpack2
235 Set(currentLevel, "CoarseNumZLayers", NumZDir);
236 Set(currentLevel, "LineDetection_Layers", TLayerId);
237 Set(currentLevel, "LineDetection_VertLineIds", TVertLineId);
238 } else {
239 Teuchos::ArrayRCP<LO> TLayerId = Teuchos::arcp<LO>(0);
240 Teuchos::ArrayRCP<LO> TVertLineId = Teuchos::arcp<LO>(0);
241 Teuchos::ArrayRCP<LO> TVertLineIdSmoo= Teuchos::arcp<LO>(0);
242
243 // store output data on current level
244 // The line detection data is used by the SemiCoarsenPFactory and the line smoothers in Ifpack/Ifpack2
245 Set(currentLevel, "CoarseNumZLayers", NumZDir);
246 Set(currentLevel, "LineDetection_Layers", TLayerId);
247 Set(currentLevel, "LineDetection_VertLineIds", TVertLineId);
248 }
249
250 // automatically switch to vertical mode on the coarser levels
251 if(Zorientation_ != VERTICAL)
252 Zorientation_ = VERTICAL;
253 }
254
255 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
256 LocalOrdinal LineDetectionFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ML_compute_line_info(LocalOrdinal LayerId[], LocalOrdinal VertLineId[], LocalOrdinal Ndof, LocalOrdinal DofsPerNode, LocalOrdinal MeshNumbering, LocalOrdinal NumNodesPerVertLine, typename Teuchos::ScalarTraits<Scalar>::coordinateType *xvals, typename Teuchos::ScalarTraits<Scalar>::coordinateType *yvals, typename Teuchos::ScalarTraits<Scalar>::coordinateType *zvals, const Teuchos::Comm<int>& /* comm */) const {
257
258 LO Nnodes, NVertLines, MyNode;
259 LO NumCoords, next; //, subindex, subnext;
260 coordinate_type xfirst, yfirst;
261 coordinate_type *xtemp, *ytemp, *ztemp;
262 LO *OrigLoc;
263 LO i,j,count;
264 LO RetVal;
265
266 RetVal = 0;
267 if ((MeshNumbering != VERTICAL) && (MeshNumbering != HORIZONTAL)) {
268 if ( (xvals == NULL) || (yvals == NULL) || (zvals == NULL)) RetVal = -1;
269 }
270 else {
271 if (NumNodesPerVertLine == -1) RetVal = -4;
272 if ( ((Ndof/DofsPerNode)%NumNodesPerVertLine) != 0) RetVal = -3;
273 }
274 if ( (Ndof%DofsPerNode) != 0) RetVal = -2;
275
276 TEUCHOS_TEST_FOR_EXCEPTION(RetVal == -1, Exceptions::RuntimeError, "Not semicoarsening as no mesh numbering information or coordinates are given\n");
277 TEUCHOS_TEST_FOR_EXCEPTION(RetVal == -4, Exceptions::RuntimeError, "Not semicoarsening as the number of z nodes is not given.\n");
278 TEUCHOS_TEST_FOR_EXCEPTION(RetVal == -3, Exceptions::RuntimeError, "Not semicoarsening as the total number of nodes is not evenly divisible by the number of z direction nodes .\n");
279 TEUCHOS_TEST_FOR_EXCEPTION(RetVal == -2, Exceptions::RuntimeError, "Not semicoarsening as something is off with the number of degrees-of-freedom per node.\n");
280
281 Nnodes = Ndof/DofsPerNode;
282 for (MyNode = 0; MyNode < Nnodes; MyNode++) VertLineId[MyNode] = -1;
283 for (MyNode = 0; MyNode < Nnodes; MyNode++) LayerId[MyNode] = -1;
284
285 if (MeshNumbering == VERTICAL) {
286 for (MyNode = 0; MyNode < Nnodes; MyNode++) {
287 LayerId[MyNode]= MyNode%NumNodesPerVertLine;
288 VertLineId[MyNode]= (MyNode- LayerId[MyNode])/NumNodesPerVertLine;
289 }
290 }
291 else if (MeshNumbering == HORIZONTAL) {
292 NVertLines = Nnodes/NumNodesPerVertLine;
293 for (MyNode = 0; MyNode < Nnodes; MyNode++) {
294 VertLineId[MyNode] = MyNode%NVertLines;
295 LayerId[MyNode] = (MyNode- VertLineId[MyNode])/NVertLines;
296 }
297 }
298 else {
299 // coordinates mode: we distinguish between vertical line numbering for semi-coarsening and line smoothing
300 NumCoords = Ndof/DofsPerNode;
301
302 // reserve temporary memory
303 Teuchos::ArrayRCP<LO> TOrigLoc= Teuchos::arcp<LO>(NumCoords); OrigLoc= TOrigLoc.getRawPtr();
304 Teuchos::ArrayRCP<coordinate_type> Txtemp = Teuchos::arcp<coordinate_type>(NumCoords); xtemp = Txtemp.getRawPtr();
305 Teuchos::ArrayRCP<coordinate_type> Tytemp = Teuchos::arcp<coordinate_type>(NumCoords); ytemp = Tytemp.getRawPtr();
306 Teuchos::ArrayRCP<coordinate_type> Tztemp = Teuchos::arcp<coordinate_type>(NumCoords); ztemp = Tztemp.getRawPtr();
307
308 // build vertical line info for semi-coarsening
309
310 // sort coordinates in {x,y,z}vals (returned in {x,y,z}temp) so that we can order things according to lines
311 // switch x and y coordinates for semi-coarsening...
312 sort_coordinates(NumCoords, OrigLoc, xvals, yvals, zvals, xtemp, ytemp, ztemp, /*true*/ true);
313
314 LO NumBlocks = 0;
315 LO index = 0;
316
317 while ( index < NumCoords ) {
318 xfirst = xtemp[index]; yfirst = ytemp[index];
319 next = index+1;
320 while ( (next != NumCoords) && (xtemp[next] == xfirst) &&
321 (ytemp[next] == yfirst))
322 next++;
323 if (NumBlocks == 0) {
324 NumNodesPerVertLine = next-index;
325 }
326 // The number of vertical lines must be the same on all processors
327 // TAW: Sep 14, 2015: or zero as we allow for empty processors.
328 //TEUCHOS_TEST_FOR_EXCEPTION(next-index != NumNodesPerVertLine,Exceptions::RuntimeError, "Error code only works for constant block size now!!!\n");
329 count = 0;
330 for (j= index; j < next; j++) {
331 VertLineId[OrigLoc[j]] = NumBlocks;
332 LayerId[OrigLoc[j]] = count++;
333 }
334 NumBlocks++;
335 index = next;
336 }
337 }
338
339 /* check that everyone was assigned */
340
341 for (i = 0; i < Nnodes; i++) {
342 if (VertLineId[i] == -1) {
343 GetOStream(Warnings1) << "Warning: did not assign " << i << " to a vertical line?????\n" << std::endl;
344 }
345 if (LayerId[i] == -1) {
346 GetOStream(Warnings1) << "Warning: did not assign " << i << " to a Layer?????\n" << std::endl;
347 }
348 }
349
350 // TAW: Sep 14 2015: relax plausibility checks as we allow for empty processors
351 //MueLu_maxAll(&comm, NumNodesPerVertLine, i);
352 //if (NumNodesPerVertLine == -1) NumNodesPerVertLine = i;
353 //TEUCHOS_TEST_FOR_EXCEPTION(NumNodesPerVertLine != i,Exceptions::RuntimeError, "Different processors have different z direction line lengths?\n");
354
355 return NumNodesPerVertLine;
356 }
357
358 /* Private member function to sort coordinates in arrays. This is an expert routine. Do not use or change.*/
359 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
361 typename Teuchos::ScalarTraits<Scalar>::coordinateType* xvals,
362 typename Teuchos::ScalarTraits<Scalar>::coordinateType* yvals,
363 typename Teuchos::ScalarTraits<Scalar>::coordinateType* zvals,
364 typename Teuchos::ScalarTraits<Scalar>::coordinateType* xtemp,
365 typename Teuchos::ScalarTraits<Scalar>::coordinateType* ytemp,
366 typename Teuchos::ScalarTraits<Scalar>::coordinateType* ztemp,
367 bool flipXY) const {
368
369 if( flipXY == false ) { // for line-smoothing
370 for (LO i = 0; i < numCoords; i++) xtemp[i]= xvals[i];
371 } else { // for semi-coarsening
372 for (LO i = 0; i < numCoords; i++) xtemp[i]= yvals[i];
373 }
374 for (LO i = 0; i < numCoords; i++) OrigLoc[i]= i;
375
376 ML_az_dsort2(xtemp,numCoords,OrigLoc);
377 if( flipXY == false ) { // for line-smoothing
378 for (LO i = 0; i < numCoords; i++) ytemp[i]= yvals[OrigLoc[i]];
379 } else {
380 for (LO i = 0; i < numCoords; i++) ytemp[i]= xvals[OrigLoc[i]];
381 }
382
383 LO index = 0;
384
385 while ( index < numCoords ) {
386 coordinate_type xfirst = xtemp[index];
387 LO next = index+1;
388 while ( (next != numCoords) && (xtemp[next] == xfirst))
389 next++;
390 ML_az_dsort2(&(ytemp[index]),next-index,&(OrigLoc[index]));
391 for (LO i = index; i < next; i++) ztemp[i]= zvals[OrigLoc[i]];
392 /* One final sort so that the ztemps are in order */
393 LO subindex = index;
394 while (subindex != next) {
395 coordinate_type yfirst = ytemp[subindex];
396 LO subnext = subindex+1;
397 while ( (subnext != next) && (ytemp[subnext] == yfirst)) subnext++;
398 ML_az_dsort2(&(ztemp[subindex]),subnext-subindex,&(OrigLoc[subindex]));
399 subindex = subnext;
400 }
401 index = next;
402 }
403
404 }
405
406 /* Sort coordinates and additional array accordingly (if provided). This is an expert routine borrowed from ML. Do not change.*/
407 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
408 void LineDetectionFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ML_az_dsort2(typename Teuchos::ScalarTraits<Scalar>::coordinateType dlist[], LocalOrdinal N, LocalOrdinal list2[]) const {
409 LO l, r, j, i, flag;
410 LO RR2;
411 coordinate_type dRR, dK;
412
413 // note: we use that routine for sorting coordinates only. No complex coordinates are assumed...
414 typedef Teuchos::ScalarTraits<SC> STS;
415
416 if (N <= 1) return;
417
418 l = N / 2 + 1;
419 r = N - 1;
420 l = l - 1;
421 dRR = dlist[l - 1];
422 dK = dlist[l - 1];
423
424 if (list2 != NULL) {
425 RR2 = list2[l - 1];
426 while (r != 0) {
427 j = l;
428 flag = 1;
429
430 while (flag == 1) {
431 i = j;
432 j = j + j;
433
434 if (j > r + 1)
435 flag = 0;
436 else {
437 if (j < r + 1)
438 if (STS::real(dlist[j]) > STS::real(dlist[j - 1])) j = j + 1;
439
440 if (STS::real(dlist[j - 1]) > STS::real(dK)) {
441 dlist[ i - 1] = dlist[ j - 1];
442 list2[i - 1] = list2[j - 1];
443 }
444 else {
445 flag = 0;
446 }
447 }
448 }
449 dlist[ i - 1] = dRR;
450 list2[i - 1] = RR2;
451
452 if (l == 1) {
453 dRR = dlist [r];
454 RR2 = list2[r];
455 dK = dlist[r];
456 dlist[r ] = dlist[0];
457 list2[r] = list2[0];
458 r = r - 1;
459 }
460 else {
461 l = l - 1;
462 dRR = dlist[ l - 1];
463 RR2 = list2[l - 1];
464 dK = dlist[l - 1];
465 }
466 }
467 dlist[ 0] = dRR;
468 list2[0] = RR2;
469 }
470 else {
471 while (r != 0) {
472 j = l;
473 flag = 1;
474 while (flag == 1) {
475 i = j;
476 j = j + j;
477 if (j > r + 1)
478 flag = 0;
479 else {
480 if (j < r + 1)
481 if (STS::real(dlist[j]) > STS::real(dlist[j - 1])) j = j + 1;
482 if (STS::real(dlist[j - 1]) > STS::real(dK)) {
483 dlist[ i - 1] = dlist[ j - 1];
484 }
485 else {
486 flag = 0;
487 }
488 }
489 }
490 dlist[ i - 1] = dRR;
491 if (l == 1) {
492 dRR = dlist [r];
493 dK = dlist[r];
494 dlist[r ] = dlist[0];
495 r = r - 1;
496 }
497 else {
498 l = l - 1;
499 dRR = dlist[ l - 1];
500 dK = dlist[l - 1];
501 }
502 }
503 dlist[ 0] = dRR;
504 }
505
506 }
507} //namespace MueLu
508
509#endif // MUELU_LINEDETECTIONFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
MueLu::DefaultLocalOrdinal LocalOrdinal
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.
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()
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void sort_coordinates(LO numCoords, LO *OrigLoc, coordinate_type *xvals, coordinate_type *yvals, coordinate_type *zvals, coordinate_type *xtemp, coordinate_type *ytemp, coordinate_type *ztemp, bool flipXY=false) const
void ML_az_dsort2(coordinate_type dlist[], LO N, LO list2[]) const
void DeclareInput(Level &currentLevel) const
Input.
typename Teuchos::ScalarTraits< SC >::coordinateType coordinate_type
void Build(Level &currentLevel) const
Build method.
LO ML_compute_line_info(LO LayerId[], LO VertLineId[], LO Ndof, LO DofsPerNode, LO MeshNumbering, LO NumNodesPerVertLine, coordinate_type *xvals, coordinate_type *yvals, coordinate_type *zvals, const Teuchos::Comm< int > &comm) const
static const NoFactory * get()
Namespace for MueLu classes and methods.
@ UserData
User data are always kept. This flag is set automatically when Level::Set("data", data) is used....
@ Runtime1
Description of what is happening (more verbose)
@ Warnings1
Additional warnings.