43 #include "PanzerAdaptersSTK_config.hpp" 48 #include "Teuchos_AbstractFactoryStd.hpp" 50 #include "Stratimikos_DefaultLinearSolverBuilder.hpp" 52 #include "Epetra_MpiComm.h" 53 #include "Epetra_Vector.h" 58 #include "Tpetra_Map.hpp" 59 #include "Tpetra_MultiVector.hpp" 61 #ifdef PANZER_HAVE_TEKO 62 #include "Teko_StratimikosFactory.hpp" 65 #ifdef PANZER_HAVE_MUELU 66 #include <Thyra_MueLuPreconditionerFactory.hpp> 67 #include <Thyra_MueLuRefMaxwellPreconditionerFactory.hpp> 68 #include "Stratimikos_MueLuHelpers.hpp" 70 #include "Xpetra_MapFactory.hpp" 71 #include "Xpetra_MultiVectorFactory.hpp" 74 #ifdef PANZER_HAVE_IFPACK2 75 #include <Thyra_Ifpack2PreconditionerFactory.hpp> 76 #include "Tpetra_CrsMatrix_decl.hpp" 86 std::vector<std::string> elementBlocks;
90 std::set<int> runningFields;
93 runningFields.insert(fields.begin(),fields.end());
97 for(std::size_t i=1;i<elementBlocks.size();i++) {
100 std::set<int> currentFields(runningFields);
101 runningFields.clear();
102 std::set_intersection(fields.begin(),fields.end(),
103 currentFields.begin(),currentFields.end(),
104 std::inserter(runningFields,runningFields.begin()));
107 if(runningFields.size()<1)
116 const std::string & fieldName,
117 std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > & fieldPatterns)
119 std::vector<std::string> elementBlocks;
122 for(std::size_t e=0;e<elementBlocks.size();e++) {
123 std::string blockId = elementBlocks[e];
126 fieldPatterns[blockId] =
127 Teuchos::rcp_dynamic_cast<const panzer::Intrepid2FieldPattern>(globalIndexer.
getFieldPattern(blockId,fieldName),
true);
133 const std::string & fieldName,
134 std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > & fieldPatterns)
137 using Teuchos::ptrFromRef;
138 using Teuchos::ptr_dynamic_cast;
143 Ptr<const DOFManager> dofManager = ptr_dynamic_cast<
const DOFManager>(ptrFromRef(globalIndexer));
145 if(dofManager!=Teuchos::null) {
146 fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
153 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
155 const Teuchos::RCP<const panzer::GlobalIndexer> & globalIndexer,
156 const Teuchos::RCP<panzer_stk::STKConnManager> & stkConn_manager,
159 const Teuchos::RCP<Teuchos::ParameterList> & strat_params,
160 #ifdef PANZER_HAVE_TEKO
161 const Teuchos::RCP<Teko::RequestHandler> & reqHandler,
163 bool writeCoordinates,
165 const Teuchos::RCP<const panzer::GlobalIndexer> & auxGlobalIndexer,
171 using Teuchos::rcp_dynamic_cast;
173 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
179 #ifdef PANZER_HAVE_MUELU 181 Stratimikos::enableMueLu<int,panzer::GlobalOrdinal,panzer::TpetraNodeType>(linearSolverBuilder,
"MueLu");
182 Stratimikos::enableMueLuRefMaxwell<int,panzer::GlobalOrdinal,panzer::TpetraNodeType>(linearSolverBuilder,
"MueLuRefMaxwell");
183 #ifndef PANZER_HIDE_DEPRECATED_CODE 185 Stratimikos::enableMueLu<int,panzer::GlobalOrdinal,panzer::TpetraNodeType>(linearSolverBuilder,
"MueLu-Tpetra");
186 Stratimikos::enableMueLuRefMaxwell<int,panzer::GlobalOrdinal,panzer::TpetraNodeType>(linearSolverBuilder,
"MueLuRefMaxwell-Tpetra");
190 #ifdef PANZER_HAVE_IFPACK2 193 typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<double, int, panzer::GlobalOrdinal,panzer::TpetraNodeType> > Impl;
195 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(),
"Ifpack2");
200 #ifdef PANZER_HAVE_TEKO 201 RCP<Teko::RequestHandler> reqHandler_local = reqHandler;
203 if(!blockedAssembly) {
205 std::string fieldName;
210 if(reqHandler_local==Teuchos::null)
211 reqHandler_local = rcp(
new Teko::RequestHandler);
214 if(determineCoordinateField(*globalIndexer,fieldName)) {
215 std::map<std::string,RCP<const panzer::Intrepid2FieldPattern> > fieldPatterns;
216 fillFieldPatternMap(*globalIndexer,fieldName,fieldPatterns);
218 RCP<panzer_stk::ParameterListCallback> callback = rcp(
new 219 panzer_stk::ParameterListCallback(fieldName,fieldPatterns,stkConn_manager,
220 rcp_dynamic_cast<const panzer::GlobalIndexer>(globalIndexer)));
221 reqHandler_local->addRequestCallback(callback);
224 if(strat_params->sublist(
"Preconditioner Types").isSublist(
"ML")) {
278 if(writeCoordinates) {
280 callback->preRequest(Teko::RequestMesg(rcp(
new Teuchos::ParameterList())));
283 const std::vector<double> & xcoords = callback->getXCoordsVector();
284 const std::vector<double> & ycoords = callback->getYCoordsVector();
285 const std::vector<double> & zcoords = callback->getZCoordsVector();
292 RCP<Epetra_Vector> vec;
307 TEUCHOS_ASSERT(
false);
311 #ifdef PANZER_HAVE_MUELU 312 if(rcp_dynamic_cast<const panzer::GlobalIndexer>(globalIndexer)!=Teuchos::null
314 if(!writeCoordinates)
315 callback->preRequest(Teko::RequestMesg(rcp(
new Teuchos::ParameterList())));
317 typedef Tpetra::Map<int,panzer::GlobalOrdinal,panzer::TpetraNodeType> Map;
318 typedef Tpetra::MultiVector<double,int,panzer::GlobalOrdinal,panzer::TpetraNodeType> MV;
321 unsigned dim = Teuchos::as<unsigned>(spatialDim);
323 for(
unsigned d=0;d<dim;d++) {
324 const std::vector<double> & coord = callback->getCoordsVector(d);
327 if(coords==Teuchos::null) {
328 if(globalIndexer->getNumFields()==1) {
329 RCP<const panzer::GlobalIndexer> ugi
331 std::vector<panzer::GlobalOrdinal> ownedIndices;
333 RCP<const Map> coords_map = rcp(
new Map(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),ownedIndices,0,mpi_comm));
334 coords = rcp(
new MV(coords_map,dim));
337 RCP<const Map> coords_map = rcp(
new Map(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),coord.size(),0,mpi_comm));
338 coords = rcp(
new MV(coords_map,dim));
343 TEUCHOS_ASSERT(coords->getLocalLength()==coord.size());
346 Teuchos::ArrayRCP<double> dest = coords->getDataNonConst(d);
347 for(std::size_t i=0;i<coord.size();i++)
352 Teuchos::ParameterList & muelu_params = strat_params->sublist(
"Preconditioner Types").sublist(
"MueLu");
353 muelu_params.set<RCP<MV> >(
"Coordinates",coords);
359 Teko::addTekoToStratimikosBuilder(linearSolverBuilder,reqHandler_local);
363 if(reqHandler_local==Teuchos::null)
364 reqHandler_local = rcp(
new Teko::RequestHandler);
366 std::string fieldName;
367 if(determineCoordinateField(*globalIndexer,fieldName)) {
368 RCP<const panzer::BlockedDOFManager> blkDofs =
370 RCP<const panzer::BlockedDOFManager> auxBlkDofs =
372 RCP<panzer_stk::ParameterListCallbackBlocked> callback =
373 rcp(
new panzer_stk::ParameterListCallbackBlocked(stkConn_manager,blkDofs,auxBlkDofs));
374 reqHandler_local->addRequestCallback(callback);
377 Teko::addTekoToStratimikosBuilder(linearSolverBuilder,reqHandler_local);
379 if(writeCoordinates) {
380 RCP<const panzer::BlockedDOFManager> blkDofs =
384 const std::vector<RCP<panzer::GlobalIndexer>> & dofVec
386 for(std::size_t i=0;i<dofVec.size();i++) {
389 TEUCHOS_ASSERT(determineCoordinateField(*dofVec[i],fieldName));
391 std::map<std::string,RCP<const panzer::Intrepid2FieldPattern> > fieldPatterns;
392 fillFieldPatternMap(*dofVec[i],fieldName,fieldPatterns);
393 panzer_stk::ParameterListCallback plCall(fieldName,fieldPatterns,stkConn_manager,dofVec[i]);
394 plCall.buildArrayToVector();
395 plCall.buildCoordinates();
398 const std::vector<double> & xcoords = plCall.getXCoordsVector();
399 const std::vector<double> & ycoords = plCall.getYCoordsVector();
400 const std::vector<double> & zcoords = plCall.getZCoordsVector();
407 RCP<Epetra_Vector> vec;
422 TEUCHOS_ASSERT(
false);
426 #ifdef PANZER_HAVE_MUELU 429 typedef Xpetra::Map<int,panzer::GlobalOrdinal> Map;
430 typedef Xpetra::MultiVector<double,int,panzer::GlobalOrdinal> MV;
433 RCP<const Map> coords_map = Xpetra::MapFactory<int,panzer::GlobalOrdinal>::Build(Xpetra::UseEpetra,
434 Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),
441 unsigned dim = Teuchos::as<unsigned>(spatialDim);
443 RCP<MV> coords = Xpetra::MultiVectorFactory<double,int,panzer::GlobalOrdinal>::Build(coords_map,spatialDim);
445 for(
unsigned d=0;d<dim;d++) {
447 TEUCHOS_ASSERT(coords->getLocalLength()==xcoords.size());
450 Teuchos::ArrayRCP<double> dest = coords->getDataNonConst(d);
451 for(std::size_t j=0;j<coords->getLocalLength();++j) {
452 if (d == 0) dest[j] = xcoords[j];
453 if (d == 1) dest[j] = ycoords[j];
454 if (d == 2) dest[j] = zcoords[j];
460 Teuchos::ParameterList & muelu_params = strat_params->sublist(
"Preconditioner Types").sublist(
"MueLu");
461 muelu_params.set<RCP<MV> >(
"Coordinates",coords);
476 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
477 "Topology writing is no longer implemented. It needs to be reimplemented for the " 478 "default DOFManager (estimate 2 days with testing)");
483 linearSolverBuilder.setParameterList(strat_params);
484 RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory = createLinearSolveStrategy(linearSolverBuilder);
490 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
492 const Teuchos::RCP<const panzer::GlobalIndexer> & globalIndexer,
493 const Teuchos::RCP<panzer::ConnManager> & conn_manager,
496 const Teuchos::RCP<Teuchos::ParameterList> & strat_params,
497 #ifdef PANZER_HAVE_TEKO
498 const Teuchos::RCP<Teko::RequestHandler> & reqHandler,
500 bool writeCoordinates,
502 const Teuchos::RCP<const panzer::GlobalIndexer> & auxGlobalIndexer,
506 #ifdef PANZER_HAVE_TEKO 507 Teuchos::RCP<Teko::RequestHandler> reqHandler_local = reqHandler;
508 if(reqHandler_local==Teuchos::null)
509 reqHandler_local = Teuchos::rcp(
new Teko::RequestHandler);
514 return buildLOWSFactory(blockedAssembly,globalIndexer,stk_conn_manager,spatialDim,mpi_comm,strat_params,
515 #ifdef PANZER_HAVE_TEKO
525 TEUCHOS_ASSERT(
false);
526 return Teuchos::null;
virtual const std::string & getFieldString(int num) const =0
Reverse lookup of the field string from a field number.
const std::vector< Teuchos::RCP< GlobalIndexer > > & getFieldDOFManagers() const
virtual void getElementBlockIds(std::vector< std::string > &elementBlockIds) const =0
void getElementBlockIds(std::vector< std::string > &elementBlockIds) const
virtual void getOwnedIndices(std::vector< panzer::GlobalOrdinal > &indices) const =0
Get the set of indices owned by this processor.
Teuchos::RCP< Thyra::LinearOpWithSolveFactoryBase< double > > buildLOWSFactory(bool blockedAssembly, const Teuchos::RCP< const panzer::GlobalIndexer > &globalIndexer, const Teuchos::RCP< panzer_stk::STKConnManager > &stkConn_manager, int spatialDim, const Teuchos::RCP< const Teuchos::MpiComm< int > > &mpi_comm, const Teuchos::RCP< Teuchos::ParameterList > &strat_params, bool writeCoordinates, bool writeTopo, const Teuchos::RCP< const panzer::GlobalIndexer > &auxGlobalIndexer, bool useCoordinates)
int VectorToMatrixMarketFile(const char *filename, const Epetra_Vector &A, const char *matrixName=0, const char *matrixDescription=0, bool writeHeader=true)
Teuchos::RCP< const FieldPattern > getFieldPattern(const std::string &name) const
Find a field pattern stored for a particular block and field number. This will retrive the pattern ad...
virtual const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const =0
bool fieldInBlock(const std::string &field, const std::string &block) const