Albany: a Trilinos-based PDE code: Albany_BCUtils

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 #include "Albany_BCUtils.hpp"
00007 
00008 // Dirichlet specialization
00009 
00010 template<>
00011 Teuchos::RCP<PHX::FieldManager<PHAL::AlbanyTraits> >
00012 Albany::BCUtils<Albany::DirichletTraits>::constructBCEvaluators(
00013   const std::vector<std::string>& nodeSetIDs,
00014   const std::vector<std::string>& bcNames,
00015   Teuchos::RCP<Teuchos::ParameterList> params,
00016   Teuchos::RCP<ParamLib> paramLib,
00017   int numEqn) {
00018 
00019   using Teuchos::RCP;
00020   using Teuchos::rcp;
00021   using Teuchos::ParameterList;
00022   using PHX::DataLayout;
00023   using PHX::MDALayout;
00024   using std::vector;
00025   using std::string;
00026 
00027   using PHAL::AlbanyTraits;
00028 
00029   if(!haveBCSpecified(params)) { // If the BC sublist is not in the input file,
00030     // but we are inside this function, this means that
00031     // node sets are contained in the Exodus file but are not defined in the problem statement.This is OK, we
00032     // just don't do anything
00033 
00034     return Teuchos::null;
00035 
00036   }
00037 
00038   Teuchos::ParameterList BCparams = params->sublist(traits_type::bcParamsPl);
00039   BCparams.validateParameters(*(traits_type::getValidBCParameters(nodeSetIDs, bcNames)), 0);
00040 
00041   std::map<std::string, RCP<ParameterList> > evaluators_to_build;
00042   RCP<DataLayout> dummy = rcp(new MDALayout<Dummy>(0));
00043   vector<std::string> bcs;
00044 
00045   // Check for all possible standard BCs (every dof on every nodeset) to see which is set
00046   for(std::size_t i = 0; i < nodeSetIDs.size(); i++) {
00047     for(std::size_t j = 0; j < bcNames.size(); j++) {
00048       std::string ss = traits_type::constructBCName(nodeSetIDs[i], bcNames[j]);
00049 
00050       if(BCparams.isParameter(ss)) {
00051         RCP<ParameterList> p = rcp(new ParameterList);
00052         p->set<int>("Type", traits_type::type);
00053 
00054         p->set< RCP<DataLayout> >("Data Layout", dummy);
00055         p->set< std::string > ("Dirichlet Name", ss);
00056         p->set< RealType >("Dirichlet Value", BCparams.get<double>(ss));
00057         p->set< std::string > ("Node Set ID", nodeSetIDs[i]);
00058         // p->set< int >     ("Number of Equations", dirichletNames.size());
00059         p->set< int > ("Equation Offset", j);
00060 
00061         p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00062 
00063         std::stringstream ess;
00064         ess << "Evaluator for " << ss;
00065         evaluators_to_build[ess.str()] = p;
00066 
00067         bcs.push_back(ss);
00068       }
00069     }
00070   }
00071 
00075   for(std::size_t i = 0; i < nodeSetIDs.size(); i++) {
00076     std::string ss = traits_type::constructBCName(nodeSetIDs[i], "CoordFunc");
00077 
00078     if(BCparams.isSublist(ss)) {
00079       // grab the sublist
00080       ParameterList& sub_list = BCparams.sublist(ss);
00081 
00082       // Directly apply the coordinate values at the boundary as a DBC (Laplace Beltrami mesh equations)
00083       if(sub_list.get<std::string>("BC Function") == "Identity") {
00084 
00085         RCP<ParameterList> p = rcp(new ParameterList);
00086         p->set<int>("Type", traits_type::typeFb);
00087 
00088         // Fill up ParameterList with things DirichletBase wants
00089         p->set< RCP<DataLayout> >("Data Layout", dummy);
00090         p->set< std::string > ("Dirichlet Name", ss);
00091         p->set< RealType > ("Dirichlet Value", 0.0);
00092         p->set< std::string > ("Node Set ID", nodeSetIDs[i]);
00093         p->set< int > ("Number of Equations", numEqn);
00094         p->set< int > ("Equation Offset", 0);
00095 
00096         p->set<RCP<ParamLib> > ("Parameter Library", paramLib);
00097         std::stringstream ess;
00098         ess << "Evaluator for " << ss;
00099         evaluators_to_build[ess.str()] = p;
00100 
00101         bcs.push_back(ss);
00102       }
00103 
00104       // Add other functional boundary conditions here. Note that Torsion could fit into this framework
00105     }
00106   }
00107 
00108 #ifdef ALBANY_LCM
00109 
00113   for(std::size_t i = 0; i < nodeSetIDs.size(); i++) {
00114     for(std::size_t j = 0; j < bcNames.size(); j++) {
00115       std::string ss = traits_type::constructTimeDepBCName(nodeSetIDs[i], bcNames[j]);
00116 
00117       if(BCparams.isSublist(ss)) {
00118         // grab the sublist
00119         ParameterList& sub_list = BCparams.sublist(ss);
00120         RCP<ParameterList> p = rcp(new ParameterList);
00121         p->set<int>("Type", traits_type::typeTd);
00122 
00123         // Extract the time values into a vector
00124         //vector<RealType> timeValues = sub_list.get<Teuchos::Array<RealType> >("Time Values").toVector();
00125         //RCP< vector<RealType> > t_ptr = Teuchos::rcpFromRef(timeValues);
00126         //p->set< RCP< vector<RealType> > >("Time Values", t_ptr);
00127         p->set< Teuchos::Array<RealType> >("Time Values", sub_list.get<Teuchos::Array<RealType> >("Time Values"));
00128 
00129         //cout << "timeValues: " << timeValues[0] << " " << timeValues[1] << std::endl;
00130 
00131         // Extract the BC values into a vector
00132         //vector<RealType> BCValues = sub_list.get<Teuchos::Array<RealType> >("BC Values").toVector();
00133         //RCP< vector<RealType> > b_ptr = Teuchos::rcpFromRef(BCValues);
00134         //p->set< RCP< vector<RealType> > >("BC Values", b_ptr);
00135         //p->set< vector<RealType> >("BC Values", BCValues);
00136         p->set< Teuchos::Array<RealType> >("BC Values", sub_list.get<Teuchos::Array<RealType> >("BC Values"));
00137         p->set< RCP<DataLayout> >("Data Layout", dummy);
00138         p->set< std::string > ("Dirichlet Name", ss);
00139         p->set< RealType >("Dirichlet Value", 0.0);
00140         p->set< int > ("Equation Offset", j);
00141         p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00142         p->set< std::string > ("Node Set ID", nodeSetIDs[i]);
00143         p->set<int>("Cubature Degree", BCparams.get("Cubature Degree", 0)); //if set to zero, the cubature degree of the side will be set to that of the element
00144 
00145         std::stringstream ess;
00146         ess << "Evaluator for " << ss;
00147         evaluators_to_build[ess.str()] = p;
00148 
00149         bcs.push_back(ss);
00150       }
00151     }
00152   }
00153 
00157   for(std::size_t i = 0; i < nodeSetIDs.size(); i++) {
00158     std::string ss = traits_type::constructBCName(nodeSetIDs[i], "twist");
00159 
00160     if(BCparams.isSublist(ss)) {
00161       // grab the sublist
00162       ParameterList& sub_list = BCparams.sublist(ss);
00163 
00164       if(sub_list.get<std::string>("BC Function") == "Torsion") {
00165         RCP<ParameterList> p = rcp(new ParameterList);
00166         p->set<int>("Type", traits_type::typeTo);
00167 
00168         p->set< RealType >("Theta Dot", sub_list.get< RealType >("Theta Dot"));
00169         p->set< RealType >("X0", sub_list.get< RealType >("X0"));
00170         p->set< RealType >("Y0", sub_list.get< RealType >("Y0"));
00171 
00172         // Fill up ParameterList with things DirichletBase wants
00173         p->set< RCP<DataLayout> >("Data Layout", dummy);
00174         p->set< std::string > ("Dirichlet Name", ss);
00175         p->set< RealType >("Dirichlet Value", 0.0);
00176         p->set< std::string > ("Node Set ID", nodeSetIDs[i]);
00177         //p->set< int >     ("Number of Equations", dirichletNames.size());
00178         p->set< int > ("Equation Offset", 0);
00179         p->set<int>("Cubature Degree", BCparams.get("Cubature Degree", 0)); //if set to zero, the cubature degree of the side will be set to that of the element
00180 
00181 
00182         p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00183         std::stringstream ess;
00184         ess << "Evaluator for " << ss;
00185         evaluators_to_build[ess.str()] = p;
00186 
00187         bcs.push_back(ss);
00188       }
00189     }
00190   }
00191 
00195   for(std::size_t i = 0; i < nodeSetIDs.size(); i++) {
00196     std::string ss = traits_type::constructBCName(nodeSetIDs[i], "K");
00197 
00198     if(BCparams.isSublist(ss)) {
00199       // grab the sublist
00200       ParameterList& sub_list = BCparams.sublist(ss);
00201 
00202       if(sub_list.get<std::string>("BC Function") == "Kfield") {
00203         RCP<ParameterList> p = rcp(new ParameterList);
00204         p->set<int>("Type", traits_type::typeKf);
00205 
00206         p->set< Teuchos::Array<RealType> >("Time Values", sub_list.get<Teuchos::Array<RealType> >("Time Values"));
00207         p->set< Teuchos::Array<RealType> >("KI Values", sub_list.get<Teuchos::Array<RealType> >("KI Values"));
00208         p->set< Teuchos::Array<RealType> >("KII Values", sub_list.get<Teuchos::Array<RealType> >("KII Values"));
00209 
00210         // // This BC needs a shear modulus and poissons ratio defined
00211         // TEUCHOS_TEST_FOR_EXCEPTION(!params->isSublist("Shear Modulus"),
00212         //               Teuchos::Exceptions::InvalidParameter,
00213         //               "This BC needs a Shear Modulus");
00214         // ParameterList& shmd_list = params->sublist("Shear Modulus");
00215         // TEUCHOS_TEST_FOR_EXCEPTION(!(shmd_list.get("Shear Modulus Type","") == "Constant"),
00216         //               Teuchos::Exceptions::InvalidParameter,
00217         //               "Invalid Shear Modulus type");
00218         // p->set< RealType >("Shear Modulus", shmd_list.get("Value", 1.0));
00219 
00220         // TEUCHOS_TEST_FOR_EXCEPTION(!params->isSublist("Poissons Ratio"),
00221         //               Teuchos::Exceptions::InvalidParameter,
00222         //               "This BC needs a Poissons Ratio");
00223         // ParameterList& pr_list = params->sublist("Poissons Ratio");
00224         // TEUCHOS_TEST_FOR_EXCEPTION(!(pr_list.get("Poissons Ratio Type","") == "Constant"),
00225         //               Teuchos::Exceptions::InvalidParameter,
00226         //               "Invalid Poissons Ratio type");
00227         // p->set< RealType >("Poissons Ratio", pr_list.get("Value", 1.0));
00228 
00229 
00230         //   p->set< Teuchos::Array<RealType> >("BC Values", sub_list.get<Teuchos::Array<RealType> >("BC Values"));
00231         //   p->set< RCP<DataLayout> >("Data Layout", dummy);
00232 
00233         // Extract BC parameters
00234         p->set< std::string >("Kfield KI Name", "Kfield KI");
00235         p->set< std::string >("Kfield KII Name", "Kfield KII");
00236         p->set< RealType >("KI Value", sub_list.get<double>("Kfield KI"));
00237         p->set< RealType >("KII Value", sub_list.get<double>("Kfield KII"));
00238         p->set< RealType >("Shear Modulus", sub_list.get<double>("Shear Modulus"));
00239         p->set< RealType >("Poissons Ratio", sub_list.get<double>("Poissons Ratio"));
00240         p->set<int>("Cubature Degree", BCparams.get("Cubature Degree", 0)); //if set to zero, the cubature degree of the side will be set to that of the element
00241 
00242 
00243         // Fill up ParameterList with things DirichletBase wants
00244         p->set< RCP<DataLayout> >("Data Layout", dummy);
00245         p->set< std::string > ("Dirichlet Name", ss);
00246         p->set< RealType >("Dirichlet Value", 0.0);
00247         p->set< std::string > ("Node Set ID", nodeSetIDs[i]);
00248         //p->set< int >     ("Number of Equations", dirichletNames.size());
00249         p->set< int > ("Equation Offset", 0);
00250 
00251         p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00252         std::stringstream ess;
00253         ess << "Evaluator for " << ss;
00254         evaluators_to_build[ess.str()] = p;
00255 
00256         bcs.push_back(ss);
00257       }
00258     }
00259   }
00260 
00261 #endif
00262 
00263   std::string allBC = "Evaluator for all Dirichlet BCs";
00264   {
00265     RCP<ParameterList> p = rcp(new ParameterList);
00266     p->set<int>("Type", traits_type::typeDa);
00267 
00268     p->set<vector<std::string>* >("DBC Names", &bcs);
00269     p->set< RCP<DataLayout> >("Data Layout", dummy);
00270     p->set<std::string>("DBC Aggregator Name", allBC);
00271     evaluators_to_build[allBC] = p;
00272   }
00273 
00274   return buildFieldManager(evaluators_to_build, allBC, dummy);
00275 }
00276 
00277 
00278 template<>
00279 Teuchos::RCP<PHX::FieldManager<PHAL::AlbanyTraits> >
00280 Albany::BCUtils<Albany::NeumannTraits>::constructBCEvaluators(
00281   const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs,
00282   const std::vector<std::string>& bcNames,
00283   const Teuchos::ArrayRCP<std::string>& dof_names,
00284   bool isVectorField,
00285   int offsetToFirstDOF,
00286   const std::vector<std::string>& conditions,
00287   const Teuchos::Array<Teuchos::Array<int> >& offsets,
00288   const Teuchos::RCP<Albany::Layouts>& dl,
00289   Teuchos::RCP<Teuchos::ParameterList> params,
00290   Teuchos::RCP<ParamLib> paramLib,
00291   const Teuchos::RCP<QCAD::MaterialDatabase>& materialDB) {
00292   using Teuchos::RCP;
00293   using Teuchos::rcp;
00294   using Teuchos::ParameterList;
00295   using PHX::DataLayout;
00296   using PHX::MDALayout;
00297   using std::vector;
00298   using std::string;
00299 
00300   using PHAL::AlbanyTraits;
00301 
00302   // Drop into the "Neumann BCs" sublist
00303   Teuchos::ParameterList BCparams = params->sublist(traits_type::bcParamsPl);
00304   BCparams.validateParameters(*(traits_type::getValidBCParameters(meshSpecs->ssNames, bcNames, conditions)), 0);
00305 
00306 
00307   std::map<std::string, RCP<ParameterList> > evaluators_to_build;
00308   vector<std::string> bcs;
00309 
00310   // Check for all possible standard BCs (every dof on every sideset) to see which is set
00311   for(std::size_t i = 0; i < meshSpecs->ssNames.size(); i++) {
00312     for(std::size_t j = 0; j < bcNames.size(); j++) {
00313       for(std::size_t k = 0; k < conditions.size(); k++) {
00314 
00315         // construct input.xml string like:
00316         // "NBC on SS sidelist_12 for DOF T set dudn"
00317         //  or
00318         // "NBC on SS sidelist_12 for DOF T set (dudx, dudy)"
00319         // or
00320         // "NBC on SS surface_1 for DOF all set P"
00321 
00322         std::string ss = traits_type::constructBCName(meshSpecs->ssNames[i], bcNames[j], conditions[k]);
00323 
00324         // Have a match of the line in input.xml
00325 
00326         if(BCparams.isParameter(ss)) {
00327 
00328           //           std::cout << "Constructing NBC: " << ss << std::endl;
00329 
00330           TEUCHOS_TEST_FOR_EXCEPTION(BCparams.isType<std::string>(ss), std::logic_error,
00331                                      "NBC array information in XML file must be of type Array(double)\n");
00332 
00333           // These are read in the Albany::Neumann constructor (PHAL_Neumann_Def.hpp)
00334 
00335           RCP<ParameterList> p = rcp(new ParameterList);
00336 
00337           p->set<int> ("Type", traits_type::type);
00338 
00339           p->set<RCP<ParamLib> > ("Parameter Library", paramLib);
00340 
00341           p->set<std::string> ("Side Set ID", meshSpecs->ssNames[i]);
00342           p->set<Teuchos::Array< int > > ("Equation Offset", offsets[j]);
00343           p->set< RCP<Albany::Layouts> > ("Layouts Struct", dl);
00344           p->set< RCP<MeshSpecsStruct> > ("Mesh Specs Struct", meshSpecs);
00345 
00346           p->set<std::string> ("Coordinate Vector Name", "Coord Vec");
00347           p->set<int>("Cubature Degree", BCparams.get("Cubature Degree", 0)); //if set to zero, the cubature degree of the side will be set to that of the element
00348 
00349           if(conditions[k] == "robin") {
00350             p->set<std::string> ("DOF Name", dof_names[j]);
00351             p->set<bool> ("Vector Field", isVectorField);
00352 
00353             if(isVectorField) {
00354               p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_vector);
00355             }
00356 
00357             else               p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_scalar);
00358           }
00359 #ifdef ALBANY_FELIX
00360           else if(conditions[k] == "basal") {
00361             std::string betaName = BCparams.get("BetaXY", "Constant");
00362             double L = BCparams.get("L", 1.0);
00363             p->set<std::string> ("BetaXY", betaName);
00364             p->set<string>("Beta Field Name", "Basal Friction");
00365             p->set<double> ("L", L);
00366             p->set<std::string> ("DOF Name", dof_names[0]);
00367             p->set<bool> ("Vector Field", isVectorField);
00368 
00369             if(isVectorField) {
00370               p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_vector);
00371             }
00372 
00373             else               p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_scalar);
00374           }
00375           else if(conditions[k] == "lateral") {
00376             std::string betaName = BCparams.get("BetaXY", "Constant");
00377             double L = BCparams.get("L", 1.0);
00378             p->set<std::string>("Thickness Field Name", "Thickness");
00379             p->set<std::string>("Elevation Field Name", "Surface Height");
00380             p->set<std::string>  ("DOF Name", dof_names[0]);
00381             p->set<bool> ("Vector Field", isVectorField);
00382             if (isVectorField) {p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_vector);}
00383             else               p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_scalar);
00384           }
00385 #endif
00386 
00387           // Pass the input file line
00388           p->set< std::string > ("Neumann Input String", ss);
00389           p->set< Teuchos::Array<double> > ("Neumann Input Value", BCparams.get<Teuchos::Array<double> >(ss));
00390           p->set< std::string > ("Neumann Input Conditions", conditions[k]);
00391 
00392           // If we are doing a Neumann internal boundary with a "scaled jump" (includes "robin" too)
00393           // The material DB database needs to be passed to the BC object
00394 
00395           if(conditions[k] == "scaled jump" || conditions[k] == "robin") {
00396 
00397             TEUCHOS_TEST_FOR_EXCEPTION(materialDB == Teuchos::null,
00398                                        Teuchos::Exceptions::InvalidParameter,
00399                                        "This BC needs a material database specified");
00400 
00401             p->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
00402 
00403 
00404           }
00405 
00406 
00407           // Inputs: X, Y at nodes, Cubature, and Basis
00408           //p->set<std::string>("Node Variable Name", "Neumann");
00409 
00410           std::stringstream ess;
00411           ess << "Evaluator for " << ss;
00412           evaluators_to_build[ess.str()] = p;
00413 
00414 
00415           bcs.push_back(ss);
00416         }
00417       }
00418     }
00419   }
00420 
00421 #ifdef ALBANY_LCM
00422 
00426   for(std::size_t i = 0; i < meshSpecs->ssNames.size(); i++) {
00427     for(std::size_t j = 0; j < bcNames.size(); j++) {
00428       for(std::size_t k = 0; k < conditions.size(); k++) {
00429 
00430         // construct input.xml string like:
00431         // "Time Dependent NBC on SS sidelist_12 for DOF T set dudn"
00432         //  or
00433         // "Time Dependent NBC on SS sidelist_12 for DOF T set (dudx, dudy)"
00434         // or
00435         // "Time Dependent NBC on SS surface_1 for DOF all set P"
00436 
00437         std::string ss = traits_type::constructTimeDepBCName(meshSpecs->ssNames[i], bcNames[j], conditions[k]);
00438 
00439         // Have a match of the line in input.xml
00440 
00441         if(BCparams.isSublist(ss)) {
00442 
00443           // grab the sublist
00444           ParameterList& sub_list = BCparams.sublist(ss);
00445 
00446           //           std::cout << "Constructing Time Dependent NBC: " << ss << std::endl;
00447 
00448           // These are read in the LCM::TimeTracBC constructor (LCM/evaluators/TimeTrac_Def.hpp)
00449 
00450           RCP<ParameterList> p = rcp(new ParameterList);
00451 
00452           p->set<int> ("Type", traits_type::typeTd);
00453 
00454           p->set< Teuchos::Array<RealType> >("Time Values",
00455                                              sub_list.get<Teuchos::Array<RealType> >("Time Values"));
00456 
00457           // Note, we use a TwoDArray here as we expect the user to specify multiple components of
00458           // the traction vector at each "time" step.
00459 
00460           p->set< Teuchos::TwoDArray<RealType> >("BC Values",
00461                                                  sub_list.get<Teuchos::TwoDArray<RealType> >("BC Values"));
00462 
00463           p->set<RCP<ParamLib> > ("Parameter Library", paramLib);
00464 
00465           p->set<std::string> ("Side Set ID", meshSpecs->ssNames[i]);
00466           p->set<Teuchos::Array< int > > ("Equation Offset", offsets[j]);
00467           p->set< RCP<Albany::Layouts> > ("Layouts Struct", dl);
00468           p->set< RCP<MeshSpecsStruct> > ("Mesh Specs Struct", meshSpecs);
00469           p->set<int>("Cubature Degree", BCparams.get("Cubature Degree", 0)); //if set to zero, the cubature degree of the side will be set to that of the element
00470 
00471 
00472           p->set<std::string> ("Coordinate Vector Name", "Coord Vec");
00473 
00474           if(conditions[k] == "robin") {
00475             p->set<std::string> ("DOF Name", dof_names[j]);
00476             p->set<bool> ("Vector Field", isVectorField);
00477 
00478             if(isVectorField) p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_vector);
00479 
00480             else               p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_scalar);
00481           }
00482 
00483           else if(conditions[k] == "basal") {
00484             p->set<std::string> ("DOF Name", dof_names[0]);
00485             p->set<bool> ("Vector Field", isVectorField);
00486 
00487             if(isVectorField) p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_vector);
00488 
00489             else               p->set< RCP<DataLayout> >("DOF Data Layout", dl->node_scalar);
00490           }
00491 
00492           // Pass the input file line
00493           p->set< std::string > ("Neumann Input String", ss);
00494           p->set< Teuchos::Array<double> > ("Neumann Input Value", Teuchos::tuple<double>(0.0, 0.0, 0.0));
00495           p->set< std::string > ("Neumann Input Conditions", conditions[k]);
00496 
00497           // If we are doing a Neumann internal boundary with a "scaled jump" (includes "robin" too)
00498           // The material DB database needs to be passed to the BC object
00499 
00500           if(conditions[k] == "scaled jump" || conditions[k] == "robin") {
00501 
00502             TEUCHOS_TEST_FOR_EXCEPTION(materialDB == Teuchos::null,
00503                                        Teuchos::Exceptions::InvalidParameter,
00504                                        "This BC needs a material database specified");
00505 
00506             p->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
00507 
00508 
00509           }
00510 
00511 
00512           std::stringstream ess;
00513           ess << "Evaluator for " << ss;
00514           evaluators_to_build[ess.str()] = p;
00515 
00516 
00517           bcs.push_back(ss);
00518         }
00519       }
00520     }
00521   }
00522 
00523 #endif
00524 
00525 
00526   // Build evaluator for Gather Coordinate Vector
00527 
00528   std::string NeuGCV = "Evaluator for Gather Coordinate Vector";
00529   {
00530     RCP<ParameterList> p = rcp(new ParameterList);
00531     p->set<int>("Type", traits_type::typeGCV);
00532 
00533     // Input: Periodic BC flag
00534     p->set<bool>("Periodic BC", false);
00535 
00536     // Output:: Coordindate Vector at vertices
00537     p->set< RCP<DataLayout> > ("Coordinate Data Layout",  dl->vertices_vector);
00538     p->set< std::string >("Coordinate Vector Name", "Coord Vec");
00539 
00540     evaluators_to_build[NeuGCV] = p;
00541   }
00542 
00543 // Build evaluator for Gather Basal Friction
00544 #ifdef ALBANY_FELIX
00545    string NeuGBF="Evaluator for Gather Basal Friction";
00546    {
00547      RCP<ParameterList> p = rcp(new ParameterList());
00548      p->set<int>("Type", traits_type::typeGBF);
00549  
00550      // for new way
00551      p->set< RCP<DataLayout> >  ("Data Layout",  dl->node_scalar);
00552      p->set< string >("Basal Friction Name", "Basal Friction");
00553 
00554      evaluators_to_build[NeuGBF] = p;
00555    }
00556 
00557    string NeuGT="Evaluator for Gather Thickness";
00558   {
00559   RCP<ParameterList> p = rcp(new ParameterList());
00560   p->set<int>("Type", traits_type::typeGT);
00561 
00562   // for new way
00563   p->set< RCP<DataLayout> >  ("Data Layout",  dl->node_scalar);
00564   p->set< string >("Thickness Name", "Thickness");
00565 
00566   evaluators_to_build[NeuGT] = p;
00567   }
00568 
00569   string NeuGSH="Evaluator for Gather Surface Height";
00570     {
00571     RCP<ParameterList> p = rcp(new ParameterList());
00572     p->set<int>("Type", traits_type::typeGSH);
00573 
00574     // for new way
00575     p->set< RCP<DataLayout> >  ("Data Layout",  dl->node_scalar);
00576     p->set< string >("Surface Height Name", "Surface Height");
00577 
00578     evaluators_to_build[NeuGSH] = p;
00579     }
00580 #endif
00581 
00582   // Build evaluator for Gather Solution
00583 
00584   std::string NeuGS = "Evaluator for Gather Solution";
00585   {
00586     RCP<ParameterList> p = rcp(new ParameterList());
00587     p->set<int>("Type", traits_type::typeGS);
00588 
00589     // for new way
00590     p->set< RCP<Albany::Layouts> >("Layouts Struct", dl);
00591 
00592     p->set< Teuchos::ArrayRCP<std::string> >("Solution Names", dof_names);
00593 
00594     p->set<bool>("Vector Field", isVectorField);
00595 
00596     if(isVectorField) p->set< RCP<DataLayout> >("Data Layout", dl->node_vector);
00597 
00598     else               p->set< RCP<DataLayout> >("Data Layout", dl->node_scalar);
00599 
00600     p->set<int>("Offset of First DOF", offsetToFirstDOF);
00601     p->set<bool>("Disable Transient", true);
00602 
00603     evaluators_to_build[NeuGS] = p;
00604   }
00605 
00606 
00607   // Build evaluator that causes the evaluation of all the NBCs
00608 
00609   std::string allBC = "Evaluator for all Neumann BCs";
00610   {
00611     RCP<ParameterList> p = rcp(new ParameterList);
00612     p->set<int>("Type", traits_type::typeNa);
00613 
00614     p->set<vector<std::string>* >("NBC Names", &bcs);
00615     p->set< RCP<DataLayout> >("Data Layout", dl->dummy);
00616     p->set<std::string>("NBC Aggregator Name", allBC);
00617     evaluators_to_build[allBC] = p;
00618   }
00619 
00620   return buildFieldManager(evaluators_to_build, allBC, dl->dummy);
00621 }
00622 
00623 template<typename BCTraits>
00624 Teuchos::RCP<PHX::FieldManager<PHAL::AlbanyTraits> >
00625 Albany::BCUtils<BCTraits>::buildFieldManager(const std::map < std::string,
00626     Teuchos::RCP<Teuchos::ParameterList> > & evaluators_to_build,
00627     std::string& allBC, Teuchos::RCP<PHX::DataLayout>& dummy) {
00628 
00629   using PHAL::AlbanyTraits;
00630 
00631   // Build Field Evaluators for each evaluation type
00632   PHX::EvaluatorFactory<AlbanyTraits, typename Albany::BCUtils<BCTraits>::traits_type::factory_type > factory;
00633 
00634   Teuchos::RCP< std::vector< Teuchos::RCP<PHX::Evaluator_TemplateManager<AlbanyTraits> > > > evaluators;
00635   evaluators = factory.buildEvaluators(evaluators_to_build);
00636 
00637   // Create a DirichletFieldManager
00638   Teuchos::RCP<PHX::FieldManager<AlbanyTraits> > fm
00639     = Teuchos::rcp(new PHX::FieldManager<AlbanyTraits>);
00640 
00641   // Register all Evaluators
00642   PHX::registerEvaluators(evaluators, *fm);
00643 
00644   PHX::Tag<AlbanyTraits::Residual::ScalarT> res_tag0(allBC, dummy);
00645   fm->requireField<AlbanyTraits::Residual>(res_tag0);
00646 
00647   PHX::Tag<AlbanyTraits::Jacobian::ScalarT> jac_tag0(allBC, dummy);
00648   fm->requireField<AlbanyTraits::Jacobian>(jac_tag0);
00649 
00650   PHX::Tag<AlbanyTraits::Tangent::ScalarT> tan_tag0(allBC, dummy);
00651   fm->requireField<AlbanyTraits::Tangent>(tan_tag0);
00652 
00653 #ifdef ALBANY_SG_MP
00654   PHX::Tag<AlbanyTraits::SGResidual::ScalarT> sgres_tag0(allBC, dummy);
00655   fm->requireField<AlbanyTraits::SGResidual>(sgres_tag0);
00656 
00657   PHX::Tag<AlbanyTraits::SGJacobian::ScalarT> sgjac_tag0(allBC, dummy);
00658   fm->requireField<AlbanyTraits::SGJacobian>(sgjac_tag0);
00659 
00660   PHX::Tag<AlbanyTraits::SGTangent::ScalarT> sgtan_tag0(allBC, dummy);
00661   fm->requireField<AlbanyTraits::SGTangent>(sgtan_tag0);
00662 
00663   PHX::Tag<AlbanyTraits::MPResidual::ScalarT> mpres_tag0(allBC, dummy);
00664   fm->requireField<AlbanyTraits::MPResidual>(mpres_tag0);
00665 
00666   PHX::Tag<AlbanyTraits::MPJacobian::ScalarT> mpjac_tag0(allBC, dummy);
00667   fm->requireField<AlbanyTraits::MPJacobian>(mpjac_tag0);
00668 
00669   PHX::Tag<AlbanyTraits::MPTangent::ScalarT> mptan_tag0(allBC, dummy);
00670   fm->requireField<AlbanyTraits::MPTangent>(mptan_tag0);
00671 #endif //ALBANY_SG_MP
00672 
00673   return fm;
00674 }
00675 
00676 // Various specializations
00677 
00678 Teuchos::RCP<const Teuchos::ParameterList>
00679 Albany::DirichletTraits::getValidBCParameters(
00680   const std::vector<std::string>& nodeSetIDs,
00681   const std::vector<std::string>& bcNames) {
00682 
00683   Teuchos::RCP<Teuchos::ParameterList> validPL =
00684     Teuchos::rcp(new Teuchos::ParameterList("Valid Dirichlet BC List"));;
00685 
00686   for(std::size_t i = 0; i < nodeSetIDs.size(); i++) {
00687     for(std::size_t j = 0; j < bcNames.size(); j++) {
00688       std::string ss = Albany::DirichletTraits::constructBCName(nodeSetIDs[i], bcNames[j]);
00689       std::string tt = Albany::DirichletTraits::constructTimeDepBCName(nodeSetIDs[i], bcNames[j]);
00690       validPL->set<double>(ss, 0.0, "Value of BC corresponding to nodeSetID and dofName");
00691       validPL->sublist(tt, false, "SubList of BC corresponding to nodeSetID and dofName");
00692     }
00693   }
00694 
00695   for(std::size_t i = 0; i < nodeSetIDs.size(); i++) {
00696     std::string ss = Albany::DirichletTraits::constructBCName(nodeSetIDs[i], "K");
00697     std::string tt = Albany::DirichletTraits::constructBCName(nodeSetIDs[i], "twist");
00698     std::string uu = Albany::DirichletTraits::constructBCName(nodeSetIDs[i], "CoordFunc");
00699     validPL->sublist(ss, false, "");
00700     validPL->sublist(tt, false, "");
00701     validPL->sublist(uu, false, "");
00702   }
00703 
00704   return validPL;
00705 
00706 }
00707 
00708 Teuchos::RCP<const Teuchos::ParameterList>
00709 Albany::NeumannTraits::getValidBCParameters(
00710   const std::vector<std::string>& sideSetIDs,
00711   const std::vector<std::string>& bcNames,
00712   const std::vector<std::string>& conditions) {
00713 
00714   Teuchos::RCP<Teuchos::ParameterList> validPL =
00715     Teuchos::rcp(new Teuchos::ParameterList("Valid Neumann BC List"));;
00716 
00717   for(std::size_t i = 0; i < sideSetIDs.size(); i++) { // loop over all side sets in the mesh
00718     for(std::size_t j = 0; j < bcNames.size(); j++) { // loop over all possible types of condition
00719       for(std::size_t k = 0; k < conditions.size(); k++) { // loop over all possible types of condition
00720 
00721         std::string ss = Albany::NeumannTraits::constructBCName(sideSetIDs[i], bcNames[j], conditions[k]);
00722         std::string tt = Albany::NeumannTraits::constructTimeDepBCName(sideSetIDs[i], bcNames[j], conditions[k]);
00723 
00724         /*
00725                 if(numDim == 2)
00726                   validPL->set<Teuchos::Array<double> >(ss, Teuchos::tuple<double>(0.0, 0.0),
00727                     "Value of BC corresponding to sideSetID and boundary condition");
00728                 else
00729                   validPL->set<Teuchos::Array<double> >(ss, Teuchos::tuple<double>(0.0, 0.0, 0.0),
00730                     "Value of BC corresponding to sideSetID and boundary condition");
00731         */
00732         Teuchos::Array<double> defaultData;
00733         validPL->set<Teuchos::Array<double> >(ss, defaultData,
00734                                               "Value of BC corresponding to sideSetID and boundary condition");
00735 
00736 
00737         validPL->sublist(tt, false, "SubList of BC corresponding to sideSetID and boundary condition");
00738       }
00739     }
00740   }
00741 
00742   validPL->set<std::string>("BetaXY", "Constant", "Function Type for Basal BC");
00743   validPL->set<int>("Cubature Degree", 3,"Cubature Degree for Neumann BC");
00744   validPL->set<double>("L", 1, "Length Scale for ISMIP-HOM Tests");
00745   return validPL;
00746 
00747 }
00748 
00749 std::string
00750 Albany::DirichletTraits::constructBCName(const std::string ns, const std::string dof) {
00751 
00752   std::stringstream ss;
00753   ss << "DBC on NS " << ns << " for DOF " << dof;
00754 
00755   return ss.str();
00756 }
00757 
00758 std::string
00759 Albany::NeumannTraits::constructBCName(const std::string ns, const std::string dof,
00760                                        const std::string condition) {
00761   std::stringstream ss;
00762   ss << "NBC on SS " << ns << " for DOF " << dof << " set " << condition;
00763   return ss.str();
00764 }
00765 
00766 std::string
00767 Albany::DirichletTraits::constructTimeDepBCName(const std::string ns, const std::string dof) {
00768   std::stringstream ss;
00769   ss << "Time Dependent " << Albany::DirichletTraits::constructBCName(ns, dof);
00770   return ss.str();
00771 }
00772 
00773 std::string
00774 Albany::NeumannTraits::constructTimeDepBCName(const std::string ns,
00775     const std::string dof, const std::string condition) {
00776   std::stringstream ss;
00777   ss << "Time Dependent " << Albany::NeumannTraits::constructBCName(ns, dof, condition);
00778   return ss.str();
00779 }
00780
Albany_BCUtils_Def.hpp
