93 FactoryMonitor m(*
this,
"Line detection (Ray style)", currentLevel);
96 RCP<CoordinateMultiVector> fineCoords;
97 ArrayRCP<coordinate_type> x, y, z;
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;
108 const ParameterList& pL = GetParameterList();
109 const std::string lineOrientation = pL.get<std::string>(
"linedetection: orientation");
113 if(lineOrientation==
"vertical")
115 else if (lineOrientation==
"horizontal")
117 else if (lineOrientation==
"coordinates")
120 TEUCHOS_TEST_FOR_EXCEPTION(
false,
Exceptions::RuntimeError,
"LineDetectionFactory: The parameter 'semicoarsen: line orientation' must be either 'vertical', 'horizontal' or 'coordinates'.");
131 GetOStream(
Runtime1) <<
"Number of layers for line detection: " << NumZDir <<
" (information from Level(0))" << std::endl;
134 NumZDir = pL.get<LO>(
"linedetection: num layers");
136 bool CoordsAvail = currentLevel.
IsAvailable(
"Coordinates");
138 if (CoordsAvail ==
true) {
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();
149 LO NumCoords = Ndofs/BlkSize;
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();
159 sort_coordinates(NumCoords, OrigLoc, xptr, yptr, zptr, xtemp, ytemp, ztemp,
true);
164 LO NumNodesPerVertLine = 0;
167 while ( index < NumCoords ) {
170 while ( (next != NumCoords) && (xtemp[next] == xfirst) &&
171 (ytemp[next] == yfirst))
173 if (NumBlocks == 0) {
174 NumNodesPerVertLine = next-index;
183 NumZDir = NumNodesPerVertLine;
184 GetOStream(
Runtime1) <<
"Number of layers for line detection: " << NumZDir <<
" (information reconstructed from provided node coordinates)" << std::endl;
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).");
189 GetOStream(
Runtime1) <<
"Number of layers for line detection: " << NumZDir <<
" (information provided by user through 'line detection: num layers')" << std::endl;
197 GetOStream(
Runtime1) <<
"Number of layers for line detection: " << NumZDir << std::endl;
199 TEUCHOS_TEST_FOR_EXCEPTION(
false,
Exceptions::RuntimeError,
"LineDetectionFactory: BuildP: No NumZLayers variable found. This cannot be.");
205 bool CoordsAvail = currentLevel.
IsAvailable(
"Coordinates");
207 if (CoordsAvail ==
false) {
211 throw Exceptions::RuntimeError(
"Coordinates not generated by previous invocation of LineDetectionFactory's BuildP() method.");
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();
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();
229 NumZDir = ML_compute_line_info(LayerId, VertLineId, Ndofs, BlkSize,
230 Zorientation_, NumZDir,xptr,yptr,zptr, *(rowMap->getComm()));
235 Set(currentLevel,
"CoarseNumZLayers", NumZDir);
236 Set(currentLevel,
"LineDetection_Layers", TLayerId);
237 Set(currentLevel,
"LineDetection_VertLineIds", TVertLineId);
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);
245 Set(currentLevel,
"CoarseNumZLayers", NumZDir);
246 Set(currentLevel,
"LineDetection_Layers", TLayerId);
247 Set(currentLevel,
"LineDetection_VertLineIds", TVertLineId);
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>& )
const {
258 LO Nnodes, NVertLines, MyNode;
268 if ( (xvals == NULL) || (yvals == NULL) || (zvals == NULL)) RetVal = -1;
271 if (NumNodesPerVertLine == -1) RetVal = -4;
272 if ( ((Ndof/DofsPerNode)%NumNodesPerVertLine) != 0) RetVal = -3;
274 if ( (Ndof%DofsPerNode) != 0) RetVal = -2;
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");
281 Nnodes = Ndof/DofsPerNode;
282 for (MyNode = 0; MyNode < Nnodes; MyNode++) VertLineId[MyNode] = -1;
283 for (MyNode = 0; MyNode < Nnodes; MyNode++) LayerId[MyNode] = -1;
286 for (MyNode = 0; MyNode < Nnodes; MyNode++) {
287 LayerId[MyNode]= MyNode%NumNodesPerVertLine;
288 VertLineId[MyNode]= (MyNode- LayerId[MyNode])/NumNodesPerVertLine;
292 NVertLines = Nnodes/NumNodesPerVertLine;
293 for (MyNode = 0; MyNode < Nnodes; MyNode++) {
294 VertLineId[MyNode] = MyNode%NVertLines;
295 LayerId[MyNode] = (MyNode- VertLineId[MyNode])/NVertLines;
300 NumCoords = Ndof/DofsPerNode;
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();
312 sort_coordinates(NumCoords, OrigLoc, xvals, yvals, zvals, xtemp, ytemp, ztemp,
true);
317 while ( index < NumCoords ) {
318 xfirst = xtemp[index]; yfirst = ytemp[index];
320 while ( (next != NumCoords) && (xtemp[next] == xfirst) &&
321 (ytemp[next] == yfirst))
323 if (NumBlocks == 0) {
324 NumNodesPerVertLine = next-index;
330 for (j= index; j < next; j++) {
331 VertLineId[OrigLoc[j]] = NumBlocks;
332 LayerId[OrigLoc[j]] = count++;
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;
345 if (LayerId[i] == -1) {
346 GetOStream(
Warnings1) <<
"Warning: did not assign " << i <<
" to a Layer?????\n" << std::endl;
355 return NumNodesPerVertLine;
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,
369 if( flipXY ==
false ) {
370 for (LO i = 0; i < numCoords; i++) xtemp[i]= xvals[i];
372 for (LO i = 0; i < numCoords; i++) xtemp[i]= yvals[i];
374 for (LO i = 0; i < numCoords; i++) OrigLoc[i]= i;
376 ML_az_dsort2(xtemp,numCoords,OrigLoc);
377 if( flipXY ==
false ) {
378 for (LO i = 0; i < numCoords; i++) ytemp[i]= yvals[OrigLoc[i]];
380 for (LO i = 0; i < numCoords; i++) ytemp[i]= xvals[OrigLoc[i]];
385 while ( index < numCoords ) {
388 while ( (next != numCoords) && (xtemp[next] == xfirst))
390 ML_az_dsort2(&(ytemp[index]),next-index,&(OrigLoc[index]));
391 for (LO i = index; i < next; i++) ztemp[i]= zvals[OrigLoc[i]];
394 while (subindex != next) {
396 LO subnext = subindex+1;
397 while ( (subnext != next) && (ytemp[subnext] == yfirst)) subnext++;
398 ML_az_dsort2(&(ztemp[subindex]),subnext-subindex,&(OrigLoc[subindex]));
414 typedef Teuchos::ScalarTraits<SC> STS;
438 if (STS::real(dlist[j]) > STS::real(dlist[j - 1])) j = j + 1;
440 if (STS::real(dlist[j - 1]) > STS::real(dK)) {
441 dlist[ i - 1] = dlist[ j - 1];
442 list2[i - 1] = list2[j - 1];
456 dlist[r ] = dlist[0];
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];
494 dlist[r ] = dlist[0];