46 #ifndef MUELU_LINEDETECTIONFACTORY_DEF_HPP 47 #define MUELU_LINEDETECTIONFACTORY_DEF_HPP 49 #include <Xpetra_Matrix.hpp> 54 #include "MueLu_FactoryManager.hpp" 61 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
63 RCP<ParameterList> validParamList = rcp(
new ParameterList());
65 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 68 #undef SET_VALID_ENTRY 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");
73 return validParamList;
76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 Input(currentLevel,
"A");
91 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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->getNodeNumElements();
105 LO Nnodes = Ndofs/BlkSize;
108 const ParameterList& pL = GetParameterList();
109 const std::string lineOrientation = pL.get<std::string>(
"linedetection: orientation");
112 if(currentLevel.GetLevelID() == 0) {
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'.");
128 if(currentLevel.GetLevelID() == 0) {
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) {
208 if (currentLevel.GetLevelID() == 0)
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);
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>& )
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;
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,
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]));
407 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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];
509 #endif // MUELU_LINEDETECTIONFACTORY_DEF_HPP #define SET_VALID_ENTRY(name)
MueLu::DefaultLocalOrdinal LocalOrdinal
void DeclareInput(Level ¤tLevel) const
Input.
Timer to be used in factories. Similar to Monitor but with additional timers.
void Build(Level ¤tLevel) const
Build method.
void ML_az_dsort2(coordinate_type dlist[], LO N, LO list2[]) const
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
typename Teuchos::ScalarTraits< SC >::coordinateType coordinate_type
Namespace for MueLu classes and methods.
static const NoFactory * get()
Class that holds all level-specific information.
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
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
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
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
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()