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" 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 "Stratimikos_MueLuHelpers.hpp" 68 #include "MatrixMarket_Tpetra.hpp" 71 #ifdef PANZER_HAVE_IFPACK2 72 #include <Thyra_Ifpack2PreconditionerFactory.hpp> 82 std::vector<std::string> elementBlocks;
86 std::set<int> runningFields;
89 runningFields.insert(fields.begin(),fields.end());
93 for(std::size_t i=1;i<elementBlocks.size();i++) {
96 std::set<int> currentFields(runningFields);
97 runningFields.clear();
98 std::set_intersection(fields.begin(),fields.end(),
99 currentFields.begin(),currentFields.end(),
100 std::inserter(runningFields,runningFields.begin()));
103 if(runningFields.size()<1)
110 #ifdef PANZER_HAVE_FEI 111 template<
typename GO>
113 fillFieldPatternMap(
const panzer::DOFManagerFEI<int,GO> & globalIndexer,
114 const std::string & fieldName,
117 std::vector<std::string> elementBlocks;
118 globalIndexer.getElementBlockIds(elementBlocks);
120 for(std::size_t e=0;e<elementBlocks.size();e++) {
121 std::string blockId = elementBlocks[e];
123 if(globalIndexer.fieldInBlock(fieldName,blockId))
124 fieldPatterns[blockId] =
125 Teuchos::rcp_dynamic_cast<const panzer::Intrepid2FieldPattern>(globalIndexer.getFieldPattern(blockId,fieldName),
true);
130 template<
typename GO>
133 const std::string & fieldName,
136 std::vector<std::string> elementBlocks;
139 for(std::size_t e=0;e<elementBlocks.size();e++) {
140 std::string blockId = elementBlocks[e];
143 fieldPatterns[blockId] =
144 Teuchos::rcp_dynamic_cast<const panzer::Intrepid2FieldPattern>(globalIndexer.
getFieldPattern(blockId,fieldName),
true);
150 const std::string & fieldName,
154 using Teuchos::ptrFromRef;
155 using Teuchos::ptr_dynamic_cast;
157 #ifdef PANZER_HAVE_FEI 158 using panzer::DOFManagerFEI;
163 Ptr<const DOFManager<int,int> > dofManager = ptr_dynamic_cast<
const DOFManager<int,int> >(ptrFromRef(globalIndexer));
165 if(dofManager!=Teuchos::null) {
166 fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
171 Ptr<const DOFManager<int,panzer::Ordinal64> > dofManager = ptr_dynamic_cast<
const DOFManager<int,panzer::Ordinal64> >(ptrFromRef(globalIndexer));
173 if(dofManager!=Teuchos::null) {
174 fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
179 #ifdef PANZER_HAVE_FEI 182 Ptr<const DOFManagerFEI<int,int> > dofManager = ptr_dynamic_cast<
const DOFManagerFEI<int,int> >(ptrFromRef(globalIndexer));
184 if(dofManager!=Teuchos::null) {
185 fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
190 Ptr<const DOFManagerFEI<int,panzer::Ordinal64> > dofManager = ptr_dynamic_cast<
const DOFManagerFEI<int,panzer::Ordinal64> >(ptrFromRef(globalIndexer));
192 if(dofManager!=Teuchos::null) {
193 fillFieldPatternMap(*dofManager,fieldName,fieldPatterns);
201 template<
typename GO>
207 const Teuchos::RCP<
const Teuchos::MpiComm<int> > & mpi_comm,
209 #ifdef PANZER_HAVE_TEKO
212 bool writeCoordinates,
216 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
222 #ifdef PANZER_HAVE_MUELU 225 Stratimikos::enableMueLu(linearSolverBuilder,
"MueLu");
226 Stratimikos::enableMueLu<int,panzer::Ordinal64,panzer::TpetraNodeType>(linearSolverBuilder,
"MueLu-Tpetra");
229 #ifdef PANZER_HAVE_IFPACK2 231 typedef Thyra::PreconditionerFactoryBase<double> Base;
232 typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<double, int, panzer::Ordinal64,panzer::TpetraNodeType> > Impl;
234 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(),
"Ifpack2");
239 #ifdef PANZER_HAVE_TEKO 242 if(!blockedAssembly) {
244 std::string fieldName;
247 if(reqHandler_local==Teuchos::null)
248 reqHandler_local =
Teuchos::rcp(
new Teko::RequestHandler);
251 if(determineCoordinateField(*globalIndexer,fieldName)) {
252 std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > fieldPatterns;
253 fillFieldPatternMap(*globalIndexer,fieldName,fieldPatterns);
256 panzer_stk::ParameterListCallback<int,GO>(fieldName,fieldPatterns,stkConn_manager,
258 reqHandler_local->addRequestCallback(callback);
260 if(writeCoordinates) {
265 const std::vector<double> & xcoords = callback->getXCoordsVector();
266 const std::vector<double> & ycoords = callback->getYCoordsVector();
267 const std::vector<double> & zcoords = callback->getZCoordsVector();
278 EpetraExt::VectorToMatrixMarketFile(
"zcoords.mm",*vec);
281 EpetraExt::VectorToMatrixMarketFile(
"ycoords.mm",*vec);
284 EpetraExt::VectorToMatrixMarketFile(
"xcoords.mm",*vec);
291 #ifdef PANZER_HAVE_MUELU 293 if(!writeCoordinates)
296 typedef Tpetra::Map<int,panzer::Ordinal64,panzer::TpetraNodeType> Map;
297 typedef Tpetra::MultiVector<double,int,panzer::Ordinal64,panzer::TpetraNodeType> MV;
300 unsigned dim = Teuchos::as<unsigned>(spatialDim);
302 for(
unsigned d=0;d<dim;d++) {
303 const std::vector<double> & coord = callback->getCoordsVector(d);
306 if(coords==Teuchos::null) {
310 std::vector<panzer::Ordinal64> ownedIndices;
311 ugi->getOwnedIndices(ownedIndices);
326 for(std::size_t i=0;i<coord.size();i++)
338 Teko::addTekoToStratimikosBuilder(linearSolverBuilder,reqHandler_local);
342 if(reqHandler_local==Teuchos::null)
343 reqHandler_local =
Teuchos::rcp(
new Teko::RequestHandler);
345 std::string fieldName;
346 if(determineCoordinateField(*globalIndexer,fieldName)) {
350 Teuchos::rcp(
new panzer_stk::ParameterListCallbackBlocked<int,GO>(stkConn_manager,blkDofs));
351 reqHandler_local->addRequestCallback(callback);
354 Teko::addTekoToStratimikosBuilder(linearSolverBuilder,reqHandler_local);
356 if(writeCoordinates) {
361 const std::vector<Teuchos::RCP<panzer::UniqueGlobalIndexer<int,GO> > > & dofVec
362 = blkDofs->getFieldDOFManagers();
363 for(std::size_t i=0;i<dofVec.size();i++) {
364 std::string fieldName;
369 std::map<std::string,Teuchos::RCP<const panzer::Intrepid2FieldPattern> > fieldPatterns;
370 fillFieldPatternMap(*dofVec[i],fieldName,fieldPatterns);
371 panzer_stk::ParameterListCallback<int,GO> plCall(fieldName,fieldPatterns,stkConn_manager,dofVec[i]);
372 plCall.buildArrayToVector();
373 plCall.buildCoordinates();
376 const std::vector<double> & xcoords = plCall.getXCoordsVector();
377 const std::vector<double> & ycoords = plCall.getYCoordsVector();
378 const std::vector<double> & zcoords = plCall.getZCoordsVector();
389 EpetraExt::VectorToMatrixMarketFile((fieldName+
"_zcoords.mm").c_str(),*vec);
392 EpetraExt::VectorToMatrixMarketFile((fieldName+
"_ycoords.mm").c_str(),*vec);
395 EpetraExt::VectorToMatrixMarketFile((fieldName+
"_xcoords.mm").c_str(),*vec);
412 linearSolverBuilder.setParameterList(strat_params);
418 template<
typename GO>
425 const std::vector<RCP<panzer::UniqueGlobalIndexer<int,GO> > > & blk_dofMngrs = blkDofs.
getFieldDOFManagers();
426 for(std::size_t b=0;b<blk_dofMngrs.size();b++) {
427 #ifdef PANZER_HAVE_FEI 428 RCP<panzer::DOFManagerFEI<int,GO> > dofMngr = Teuchos::rcp_dynamic_cast<panzer::DOFManagerFEI<int,GO> >(blk_dofMngrs[b],
true);
430 std::vector<std::string> eBlocks;
431 dofMngr->getElementBlockIds(eBlocks);
434 std::stringstream fileName;
435 fileName <<
"elements_" << b;
436 std::ofstream file(fileName.str().c_str());
439 for(std::size_t e=0;e<eBlocks.size();e++)
447 #ifdef PANZER_HAVE_FEI 448 template <
typename GO>
450 writeTopology(
const panzer::DOFManagerFEI<int,GO> & dofs,
const std::string & block,std::ostream & os)
452 std::vector<std::string> fields(dofs.getElementBlockGIDCount(block));
454 const std::set<int> & fieldIds = dofs.getFields(block);
455 for(std::set<int>::const_iterator itr=fieldIds.begin();itr!=fieldIds.end();++itr) {
456 std::string
field = dofs.getFieldString(*itr);
459 const std::vector<int> & fieldOffsets = dofs.getGIDFieldOffsets(block,*itr);
460 for(std::size_t f=0;f<fieldOffsets.size();f++)
461 fields[fieldOffsets[f]] =
field;
466 os <<
"#" << std::endl;
467 os <<
"# Element Block \"" << block <<
"\"" << std::endl;
468 os <<
"# field pattern = [ " << fields[0];
469 for(std::size_t f=1;f<fields.size();f++)
470 os <<
", " << fields[f];
471 os <<
" ]" << std::endl;
472 os <<
"#" << std::endl;
474 const std::vector<int> & elements = dofs.getElementBlock(block);
475 for(std::size_t e=0;e<elements.size();e++) {
476 std::vector<GO> gids;
477 dofs.getElementGIDs(elements[e],gids,block);
480 os <<
"[ " << gids[0];
481 for(std::size_t g=1;g<gids.size();g++)
482 os <<
", " << gids[g];
483 os <<
" ]" << std::endl;
virtual int getNumFields() const =0
void getElementBlockIds(std::vector< std::string > &elementBlockIds) const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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
PHX::MDField< const ScalarT, Cell, IP > field
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
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
#define TEUCHOS_ASSERT(assertion_test)
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)
void writeTopology(const panzer::BlockedDOFManager< int, GO > &blkDofs)