43 #ifndef __Panzer_STK_SetupLOWSFactory_impl_hpp__ 44 #define __Panzer_STK_SetupLOWSFactory_impl_hpp__ 46 #include "PanzerAdaptersSTK_config.hpp" 50 #include "Teuchos_AbstractFactoryStd.hpp" 52 #include "Stratimikos_DefaultLinearSolverBuilder.hpp" 54 #include "Epetra_MpiComm.h" 55 #include "Epetra_Vector.h" 56 #include "EpetraExt_VectorOut.h" 60 #include "Tpetra_Map.hpp" 61 #include "Tpetra_MultiVector.hpp" 63 #ifdef PANZER_HAVE_TEKO 64 #include "Teko_StratimikosFactory.hpp" 67 #ifdef PANZER_HAVE_MUELU 68 #include <Thyra_MueLuPreconditionerFactory.hpp> 69 #include <Thyra_MueLuRefMaxwellPreconditionerFactory.hpp> 70 #include "Stratimikos_MueLuHelpers.hpp" 71 #include "MatrixMarket_Tpetra.hpp" 72 #include "Xpetra_MapFactory.hpp" 73 #include "Xpetra_MultiVectorFactory.hpp" 76 #ifdef PANZER_HAVE_IFPACK2 77 #include <Thyra_Ifpack2PreconditionerFactory.hpp> 87 std::vector<std::string> elementBlocks;
91 std::set<int> runningFields;
94 runningFields.insert(fields.begin(),fields.end());
98 for(std::size_t i=1;i<elementBlocks.size();i++) {
101 std::set<int> currentFields(runningFields);
102 runningFields.clear();
103 std::set_intersection(fields.begin(),fields.end(),
104 currentFields.begin(),currentFields.end(),
105 std::inserter(runningFields,runningFields.begin()));
108 if(runningFields.size()<1)
115 template<
typename GO>
118 const std::string & fieldName,
119 std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > & fieldPatterns)
121 std::vector<std::string> elementBlocks;
124 for(std::size_t e=0;e<elementBlocks.size();e++) {
125 std::string blockId = elementBlocks[e];
128 fieldPatterns[blockId] =
129 Teuchos::rcp_dynamic_cast<const panzer::Intrepid2FieldPattern>(globalIndexer.
getFieldPattern(blockId,fieldName),
true);
135 const std::string & fieldName,
136 std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > & fieldPatterns)
139 using Teuchos::ptrFromRef;
140 using Teuchos::ptr_dynamic_cast;
145 Ptr<const DOFManager<int,int> > dofManager = ptr_dynamic_cast<
const DOFManager<int,int> >(ptrFromRef(globalIndexer));
147 if(dofManager!=Teuchos::null) {
148 fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
153 Ptr<const DOFManager<int,panzer::Ordinal64> > dofManager = ptr_dynamic_cast<
const DOFManager<int,panzer::Ordinal64> >(ptrFromRef(globalIndexer));
155 if(dofManager!=Teuchos::null) {
156 fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
163 template<
typename GO>
164 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
166 const Teuchos::RCP<const panzer::UniqueGlobalIndexerBase> & globalIndexer,
170 const Teuchos::RCP<Teuchos::ParameterList> & strat_params,
171 #ifdef PANZER_HAVE_TEKO
172 const Teuchos::RCP<Teko::RequestHandler> & reqHandler,
174 bool writeCoordinates,
176 const Teuchos::RCP<const panzer::UniqueGlobalIndexerBase> & auxGlobalIndexer,
182 using Teuchos::rcp_dynamic_cast;
184 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
190 #ifdef PANZER_HAVE_MUELU 193 Stratimikos::enableMueLu(linearSolverBuilder,
"MueLu");
194 Stratimikos::enableMueLuRefMaxwell(linearSolverBuilder,
"MueLuRefMaxwell");
195 Stratimikos::enableMueLu<int,panzer::Ordinal64,panzer::TpetraNodeType>(linearSolverBuilder,
"MueLu-Tpetra");
196 Stratimikos::enableMueLuRefMaxwell<int,panzer::Ordinal64,panzer::TpetraNodeType>(linearSolverBuilder,
"MueLuRefMaxwell-Tpetra");
199 #ifdef PANZER_HAVE_IFPACK2 201 typedef Thyra::PreconditionerFactoryBase<double> Base;
202 typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<double, int, panzer::Ordinal64,panzer::TpetraNodeType> > Impl;
204 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(),
"Ifpack2");
209 #ifdef PANZER_HAVE_TEKO 210 RCP<Teko::RequestHandler> reqHandler_local = reqHandler;
212 if(!blockedAssembly) {
214 std::string fieldName;
219 if(reqHandler_local==Teuchos::null)
220 reqHandler_local = rcp(
new Teko::RequestHandler);
223 if(determineCoordinateField(*globalIndexer,fieldName)) {
224 std::map<std::string,RCP<const panzer::Intrepid2FieldPattern> > fieldPatterns;
225 fillFieldPatternMap(*globalIndexer,fieldName,fieldPatterns);
227 RCP<panzer_stk::ParameterListCallback<int,GO> > callback = rcp(
new 228 panzer_stk::ParameterListCallback<int,GO>(fieldName,fieldPatterns,stkConn_manager,
230 reqHandler_local->addRequestCallback(callback);
233 if(strat_params->sublist(
"Preconditioner Types").isSublist(
"ML")) {
287 if(writeCoordinates) {
289 callback->preRequest(Teko::RequestMesg(rcp(
new Teuchos::ParameterList())));
292 const std::vector<double> & xcoords = callback->getXCoordsVector();
293 const std::vector<double> & ycoords = callback->getYCoordsVector();
294 const std::vector<double> & zcoords = callback->getZCoordsVector();
297 Epetra_MpiComm ep_comm(*mpi_comm->getRawMpiComm());
299 Epetra_Map map(-1,xcoords.size(),0,ep_comm);
301 RCP<Epetra_Vector> vec;
304 vec = rcp(
new Epetra_Vector(Copy,map,const_cast<double *>(&zcoords[0])));
305 EpetraExt::VectorToMatrixMarketFile(
"zcoords.mm",*vec);
308 vec = rcp(
new Epetra_Vector(Copy,map,const_cast<double *>(&ycoords[0])));
309 EpetraExt::VectorToMatrixMarketFile(
"ycoords.mm",*vec);
312 vec = rcp(
new Epetra_Vector(Copy,map,const_cast<double *>(&xcoords[0])));
313 EpetraExt::VectorToMatrixMarketFile(
"xcoords.mm",*vec);
316 TEUCHOS_ASSERT(
false);
320 #ifdef PANZER_HAVE_MUELU 323 if(!writeCoordinates)
324 callback->preRequest(Teko::RequestMesg(rcp(
new Teuchos::ParameterList())));
326 typedef Tpetra::Map<int,panzer::Ordinal64,panzer::TpetraNodeType> Map;
327 typedef Tpetra::MultiVector<double,int,panzer::Ordinal64,panzer::TpetraNodeType> MV;
330 unsigned dim = Teuchos::as<unsigned>(spatialDim);
332 for(
unsigned d=0;d<dim;d++) {
333 const std::vector<double> & coord = callback->getCoordsVector(d);
336 if(coords==Teuchos::null) {
337 if(globalIndexer->getNumFields()==1) {
338 RCP<const panzer::UniqueGlobalIndexer<int,panzer::Ordinal64> > ugi
340 std::vector<panzer::Ordinal64> ownedIndices;
342 RCP<const Map> coords_map = rcp(
new Map(Teuchos::OrdinalTraits<panzer::Ordinal64>::invalid(),ownedIndices,0,mpi_comm));
343 coords = rcp(
new MV(coords_map,dim));
346 RCP<const Map> coords_map = rcp(
new Map(Teuchos::OrdinalTraits<panzer::Ordinal64>::invalid(),coord.size(),0,mpi_comm));
347 coords = rcp(
new MV(coords_map,dim));
352 TEUCHOS_ASSERT(coords->getLocalLength()==coord.size());
355 Teuchos::ArrayRCP<double> dest = coords->getDataNonConst(d);
356 for(std::size_t i=0;i<coord.size();i++)
361 Teuchos::ParameterList & muelu_params = strat_params->sublist(
"Preconditioner Types").sublist(
"MueLu-Tpetra");
362 muelu_params.set<RCP<MV> >(
"Coordinates",coords);
368 Teko::addTekoToStratimikosBuilder(linearSolverBuilder,reqHandler_local);
372 if(reqHandler_local==Teuchos::null)
373 reqHandler_local = rcp(
new Teko::RequestHandler);
375 std::string fieldName;
376 if(determineCoordinateField(*globalIndexer,fieldName)) {
377 RCP<const panzer::BlockedDOFManager<int,GO> > blkDofs =
379 RCP<const panzer::BlockedDOFManager<int,GO> > auxBlkDofs =
381 RCP<panzer_stk::ParameterListCallbackBlocked<int,GO> > callback =
382 rcp(
new panzer_stk::ParameterListCallbackBlocked<int,GO>(stkConn_manager,blkDofs,auxBlkDofs));
383 reqHandler_local->addRequestCallback(callback);
386 Teko::addTekoToStratimikosBuilder(linearSolverBuilder,reqHandler_local);
388 if(writeCoordinates) {
389 RCP<const panzer::BlockedDOFManager<int,GO> > blkDofs =
393 const std::vector<RCP<panzer::UniqueGlobalIndexer<int,GO> > > & dofVec
395 for(std::size_t i=0;i<dofVec.size();i++) {
396 std::string fieldName;
399 TEUCHOS_ASSERT(determineCoordinateField(*dofVec[i],fieldName));
401 std::map<std::string,RCP<const panzer::Intrepid2FieldPattern> > fieldPatterns;
402 fillFieldPatternMap(*dofVec[i],fieldName,fieldPatterns);
403 panzer_stk::ParameterListCallback<int,GO> plCall(fieldName,fieldPatterns,stkConn_manager,dofVec[i]);
404 plCall.buildArrayToVector();
405 plCall.buildCoordinates();
408 const std::vector<double> & xcoords = plCall.getXCoordsVector();
409 const std::vector<double> & ycoords = plCall.getYCoordsVector();
410 const std::vector<double> & zcoords = plCall.getZCoordsVector();
413 Epetra_MpiComm ep_comm(*mpi_comm->getRawMpiComm());
415 Epetra_Map map(-1,xcoords.size(),0,ep_comm);
417 RCP<Epetra_Vector> vec;
420 vec = rcp(
new Epetra_Vector(Copy,map,const_cast<double *>(&zcoords[0])));
421 EpetraExt::VectorToMatrixMarketFile((fieldName+
"_zcoords.mm").c_str(),*vec);
424 vec = rcp(
new Epetra_Vector(Copy,map,const_cast<double *>(&ycoords[0])));
425 EpetraExt::VectorToMatrixMarketFile((fieldName+
"_ycoords.mm").c_str(),*vec);
428 vec = rcp(
new Epetra_Vector(Copy,map,const_cast<double *>(&xcoords[0])));
429 EpetraExt::VectorToMatrixMarketFile((fieldName+
"_xcoords.mm").c_str(),*vec);
432 TEUCHOS_ASSERT(
false);
436 #ifdef PANZER_HAVE_MUELU 439 typedef Xpetra::Map<int,GO> Map;
440 typedef Xpetra::MultiVector<double,int,GO> MV;
443 RCP<const Map> coords_map = Xpetra::MapFactory<int,GO>::Build(Xpetra::UseEpetra,
444 Teuchos::OrdinalTraits<GO>::invalid(),
451 unsigned dim = Teuchos::as<unsigned>(spatialDim);
453 RCP<MV> coords = Xpetra::MultiVectorFactory<double,int,GO>::Build(coords_map,spatialDim);
455 for(
unsigned d=0;d<dim;d++) {
457 TEUCHOS_ASSERT(coords->getLocalLength()==xcoords.size());
460 Teuchos::ArrayRCP<double> dest = coords->getDataNonConst(d);
461 for(std::size_t i=0;i<coords->getLocalLength();i++) {
462 if (d == 0) dest[i] = xcoords[i];
463 if (d == 1) dest[i] = ycoords[i];
464 if (d == 2) dest[i] = zcoords[i];
470 Teuchos::ParameterList & muelu_params = strat_params->sublist(
"Preconditioner Types").sublist(
"MueLu");
471 muelu_params.set<RCP<MV> >(
"Coordinates",coords);
486 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
487 "Topology writing is no longer implemented. It needs to be reimplemented for the " 488 "default DOFManager (estimate 2 days with testing)");
493 linearSolverBuilder.setParameterList(strat_params);
494 RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory = createLinearSolveStrategy(linearSolverBuilder);
Teuchos::RCP< Thyra::LinearOpWithSolveFactoryBase< double > > buildLOWSFactory(bool blockedAssembly, const Teuchos::RCP< const panzer::UniqueGlobalIndexerBase > &globalIndexer, const Teuchos::RCP< panzer::ConnManagerBase< int > > &conn_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::UniqueGlobalIndexerBase > &auxGlobalIndexer, bool useCoordinates)
void getElementBlockIds(std::vector< std::string > &elementBlockIds) const
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 void getElementBlockIds(std::vector< std::string > &elementBlockIds) const =0
virtual const std::vector< int > & getBlockFieldNumbers(const std::string &blockId) const =0
const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > & getFieldDOFManagers() const
virtual const std::string & getFieldString(int num) const =0
Reverse lookup of the field string from a field number.
bool fieldInBlock(const std::string &field, const std::string &block) const
virtual void getOwnedIndices(std::vector< GlobalOrdinalT > &indices) const =0
Get the set of indices owned by this processor.