46 #ifndef MUELU_LINEDETECTIONFACTORY_DEF_HPP 47 #define MUELU_LINEDETECTIONFACTORY_DEF_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<MultiVector> fineCoords;
97 ArrayRCP<Scalar> x, y, z;
98 Scalar *xptr = NULL, *yptr = NULL, *zptr = NULL;
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<MultiVector> > (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<SC> Txtemp = Teuchos::arcp<SC>(NumCoords);
SC* xtemp = Txtemp.getRawPtr();
154 Teuchos::ArrayRCP<SC> Tytemp = Teuchos::arcp<SC>(NumCoords);
SC* ytemp = Tytemp.getRawPtr();
155 Teuchos::ArrayRCP<SC> Tztemp = Teuchos::arcp<SC>(NumCoords);
SC* ztemp = Tztemp.getRawPtr();
159 sort_coordinates(NumCoords, OrigLoc, xptr, yptr, zptr, xtemp, ytemp, ztemp,
true);
164 LO NumNodesPerVertLine = 0;
167 while ( index < NumCoords ) {
168 SC xfirst = xtemp[index];
SC yfirst = ytemp[index];
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<MultiVector> > (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, Scalar *xvals, Scalar *yvals, Scalar *zvals,
const Teuchos::Comm<int>& comm)
const {
258 LO Nnodes, NVertLines, MyNode;
261 SC *xtemp, *ytemp, *ztemp;
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<SC> Txtemp = Teuchos::arcp<SC>(NumCoords); xtemp = Txtemp.getRawPtr();
305 Teuchos::ArrayRCP<SC> Tytemp = Teuchos::arcp<SC>(NumCoords); ytemp = Tytemp.getRawPtr();
306 Teuchos::ArrayRCP<SC> Tztemp = Teuchos::arcp<SC>(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>
360 void LineDetectionFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sort_coordinates(
LO numCoords,
LO* OrigLoc, Scalar* xvals, Scalar* yvals, Scalar* zvals, Scalar* xtemp, Scalar* ytemp, Scalar* ztemp,
bool flipXY)
const {
362 if( flipXY ==
false ) {
363 for (
LO i = 0; i < numCoords; i++) xtemp[i]= xvals[i];
365 for (
LO i = 0; i < numCoords; i++) xtemp[i]= yvals[i];
367 for (
LO i = 0; i < numCoords; i++) OrigLoc[i]= i;
369 ML_az_dsort2(xtemp,numCoords,OrigLoc);
370 if( flipXY ==
false ) {
371 for (
LO i = 0; i < numCoords; i++) ytemp[i]= yvals[OrigLoc[i]];
373 for (
LO i = 0; i < numCoords; i++) ytemp[i]= xvals[OrigLoc[i]];
378 while ( index < numCoords ) {
379 SC xfirst = xtemp[index];
381 while ( (next != numCoords) && (xtemp[next] == xfirst))
383 ML_az_dsort2(&(ytemp[index]),next-index,&(OrigLoc[index]));
384 for (
LO i = index; i < next; i++) ztemp[i]= zvals[OrigLoc[i]];
387 while (subindex != next) {
388 SC yfirst = ytemp[subindex];
389 LO subnext = subindex+1;
390 while ( (subnext != next) && (ytemp[subnext] == yfirst)) subnext++;
391 ML_az_dsort2(&(ztemp[subindex]),subnext-subindex,&(OrigLoc[subindex]));
400 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
407 typedef Teuchos::ScalarTraits<SC> STS;
431 if (STS::real(dlist[j]) > STS::real(dlist[j - 1])) j = j + 1;
433 if (STS::real(dlist[j - 1]) > STS::real(dK)) {
434 dlist[ i - 1] = dlist[ j - 1];
435 list2[i - 1] = list2[j - 1];
449 dlist[r ] = dlist[0];
474 if (STS::real(dlist[j]) > STS::real(dlist[j - 1])) j = j + 1;
475 if (STS::real(dlist[j - 1]) > STS::real(dK)) {
476 dlist[ i - 1] = dlist[ j - 1];
487 dlist[r ] = dlist[0];
502 #endif // MUELU_LINEDETECTIONFACTORY_DEF_HPP #define SET_VALID_ENTRY(name)
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.
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
Namespace for MueLu classes and methods.
void sort_coordinates(LO numCoords, LO *OrigLoc, Scalar *xvals, Scalar *yvals, Scalar *zvals, Scalar *xtemp, Scalar *ytemp, Scalar *ztemp, bool flipXY=false) const
static const NoFactory * get()
Class that holds all level-specific information.
LO ML_compute_line_info(LO LayerId[], LO VertLineId[], LO Ndof, LO DofsPerNode, LO MeshNumbering, LO NumNodesPerVertLine, SC *xvals, SC *yvals, SC *zvals, const Teuchos::Comm< int > &comm) const
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
void ML_az_dsort2(SC dlist[], LO N, LO list2[]) 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)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()