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<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 sort_coordinates(LO numCoords, LO *OrigLoc, Scalar *xvals, Scalar *yvals, Scalar *zvals, Scalar *xtemp, Scalar *ytemp, Scalar *ztemp, bool flipXY=false) const
Timer to be used in factories. Similar to Monitor but with additional timers.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
User data are always kept. This flag is set automatically when Level::Set("data", data) is used...
Namespace for MueLu classes and methods.
static const NoFactory * get()
void Build(Level ¤tLevel) const
Build method.
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
Class that holds all level-specific information.
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
void ML_az_dsort2(SC dlist[], LO N, LO list2[]) const
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()
void DeclareInput(Level ¤tLevel) const
Input.