00001
00002
00003
00004
00005
00006
00007 #ifndef FELIX_STOKESFOPROBLEM_HPP
00008 #define FELIX_STOKESFOPROBLEM_HPP
00009
00010 #include "Teuchos_RCP.hpp"
00011 #include "Teuchos_ParameterList.hpp"
00012
00013 #include "Albany_AbstractProblem.hpp"
00014
00015 #include "Phalanx.hpp"
00016 #include "PHAL_Workset.hpp"
00017 #include "PHAL_Dimension.hpp"
00018
00019
00020
00021
00022 namespace FELIX {
00023
00028 class StokesFO : public Albany::AbstractProblem {
00029 public:
00030
00032 StokesFO(const Teuchos::RCP<Teuchos::ParameterList>& params,
00033 const Teuchos::RCP<ParamLib>& paramLib,
00034 const int numDim_);
00035
00037 ~StokesFO();
00038
00040 virtual int spatialDimension() const { return numDim; }
00041
00043 virtual void buildProblem(
00044 Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> > meshSpecs,
00045 Albany::StateManager& stateMgr);
00046
00047
00048 virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00049 buildEvaluators(
00050 PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00051 const Albany::MeshSpecsStruct& meshSpecs,
00052 Albany::StateManager& stateMgr,
00053 Albany::FieldManagerChoice fmchoice,
00054 const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00055
00057 Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00058
00059 private:
00060
00062 StokesFO(const StokesFO&);
00063
00065 StokesFO& operator=(const StokesFO&);
00066
00067 public:
00068
00070 template <typename EvalT> Teuchos::RCP<const PHX::FieldTag>
00071 constructEvaluators(
00072 PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00073 const Albany::MeshSpecsStruct& meshSpecs,
00074 Albany::StateManager& stateMgr,
00075 Albany::FieldManagerChoice fmchoice,
00076 const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00077
00078 void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00079 void constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs);
00080
00081 protected:
00082 int numDim;
00083 Teuchos::RCP<Albany::Layouts> dl;
00084
00085 };
00086
00087 }
00088
00089 #include "Intrepid_FieldContainer.hpp"
00090 #include "Intrepid_DefaultCubatureFactory.hpp"
00091 #include "Shards_CellTopology.hpp"
00092
00093 #include "Albany_Utils.hpp"
00094 #include "Albany_ProblemUtils.hpp"
00095 #include "Albany_EvaluatorUtils.hpp"
00096 #include "Albany_ResponseUtilities.hpp"
00097
00098 #include "FELIX_StokesFOResid.hpp"
00099 #include "FELIX_ViscosityFO.hpp"
00100 #include "FELIX_StokesFOBodyForce.hpp"
00101 #include "PHAL_Neumann.hpp"
00102 #include "PHAL_Source.hpp"
00103
00104
00105 template <typename EvalT>
00106 Teuchos::RCP<const PHX::FieldTag>
00107 FELIX::StokesFO::constructEvaluators(
00108 PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00109 const Albany::MeshSpecsStruct& meshSpecs,
00110 Albany::StateManager& stateMgr,
00111 Albany::FieldManagerChoice fieldManagerChoice,
00112 const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00113 {
00114 using Teuchos::RCP;
00115 using Teuchos::rcp;
00116 using Teuchos::ParameterList;
00117 using PHX::DataLayout;
00118 using PHX::MDALayout;
00119 using std::vector;
00120 using std::string;
00121 using std::map;
00122 using PHAL::AlbanyTraits;
00123
00124 RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00125 intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00126 RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00127
00128 const int numNodes = intrepidBasis->getCardinality();
00129 const int worksetSize = meshSpecs.worksetSize;
00130
00131 Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00132 RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00133
00134 const int numQPts = cubature->getNumPoints();
00135 const int numVertices = cellType->getNodeCount();
00136 int vecDim = neq;
00137
00138 #ifdef OUTPUT_TO_SCREEN
00139 *out << "Field Dimensions: Workset=" << worksetSize
00140 << ", Vertices= " << numVertices
00141 << ", Nodes= " << numNodes
00142 << ", QuadPts= " << numQPts
00143 << ", Dim= " << numDim
00144 << ", vecDim= " << vecDim << std::endl;
00145 #endif
00146
00147 dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim, vecDim));
00148 Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00149 int offset=0;
00150
00151
00152 Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00153
00154
00155
00156 Teuchos::ArrayRCP<std::string> dof_names(1);
00157 Teuchos::ArrayRCP<std::string> dof_names_dot(1);
00158 Teuchos::ArrayRCP<std::string> resid_names(1);
00159 dof_names[0] = "Velocity";
00160
00161 resid_names[0] = "Stokes Residual";
00162 fm0.template registerEvaluator<EvalT>
00163 (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, offset));
00164
00165
00166 fm0.template registerEvaluator<EvalT>
00167 (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0]));
00168
00169
00170
00171
00172 fm0.template registerEvaluator<EvalT>
00173 (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0]));
00174
00175 fm0.template registerEvaluator<EvalT>
00176 (evalUtils.constructScatterResidualEvaluator(true, resid_names,offset, "Scatter Stokes"));
00177 offset += numDim;
00178
00179 fm0.template registerEvaluator<EvalT>
00180 (evalUtils.constructGatherCoordinateVectorEvaluator());
00181
00182 fm0.template registerEvaluator<EvalT>
00183 (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00184
00185 fm0.template registerEvaluator<EvalT>
00186 (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00187
00188 fm0.template registerEvaluator<EvalT>
00189 (evalUtils.constructGatherSHeightEvaluator());
00190
00191 fm0.template registerEvaluator<EvalT>
00192 (evalUtils.constructGatherTemperatureEvaluator());
00193
00194 fm0.template registerEvaluator<EvalT>
00195 (evalUtils.constructGatherFlowFactorEvaluator());
00196
00197 std::string sh = "Surface Height";
00198 fm0.template registerEvaluator<EvalT>
00199 (evalUtils.constructDOFGradInterpolationEvaluator_noDeriv(sh));
00200
00201 {
00202 RCP<ParameterList> p = rcp(new ParameterList("Stokes Resid"));
00203
00204
00205 p->set<std::string>("Weighted BF Name", "wBF");
00206 p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00207 p->set<std::string>("QP Variable Name", "Velocity");
00208 p->set<std::string>("QP Time Derivative Variable Name", "Velocity_dot");
00209 p->set<std::string>("Gradient QP Variable Name", "Velocity Gradient");
00210 p->set<std::string>("Body Force Name", "Body Force");
00211 p->set<std::string>("FELIX Viscosity QP Variable Name", "FELIX Viscosity");
00212
00213 Teuchos::ParameterList& paramList = params->sublist("Equation Set");
00214 p->set<Teuchos::ParameterList*>("Parameter List", ¶mList);
00215
00216
00217 p->set<std::string>("Residual Name", "Stokes Residual");
00218
00219 ev = rcp(new FELIX::StokesFOResid<EvalT,AlbanyTraits>(*p,dl));
00220 fm0.template registerEvaluator<EvalT>(ev);
00221 }
00222 {
00223 RCP<ParameterList> p = rcp(new ParameterList("FELIX Viscosity"));
00224
00225
00226 p->set<std::string>("Coordinate Vector Name", "Coord Vec");
00227 p->set<std::string>("Gradient QP Variable Name", "Velocity Gradient");
00228 p->set<std::string>("Temperature Name", "Temperature");
00229 p->set<std::string>("Flow Factor Name", "Flow Factor");
00230
00231 p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00232 Teuchos::ParameterList& paramList = params->sublist("FELIX Viscosity");
00233 p->set<Teuchos::ParameterList*>("Parameter List", ¶mList);
00234
00235
00236 p->set<std::string>("FELIX Viscosity QP Variable Name", "FELIX Viscosity");
00237
00238 ev = rcp(new FELIX::ViscosityFO<EvalT,AlbanyTraits>(*p,dl));
00239 fm0.template registerEvaluator<EvalT>(ev);
00240
00241 }
00242
00243 {
00244 RCP<ParameterList> p = rcp(new ParameterList("Body Force"));
00245
00246
00247 p->set<std::string>("FELIX Viscosity QP Variable Name", "FELIX Viscosity");
00248 p->set<std::string>("Coordinate Vector Name", "Coord Vec");
00249 p->set<std::string>("Surface Height Gradient Name", "Surface Height Gradient");
00250
00251 Teuchos::ParameterList& paramList = params->sublist("Body Force");
00252 p->set<Teuchos::ParameterList*>("Parameter List", ¶mList);
00253
00254
00255
00256 p->set<std::string>("Body Force Name", "Body Force");
00257
00258 ev = rcp(new FELIX::StokesFOBodyForce<EvalT,AlbanyTraits>(*p,dl));
00259 fm0.template registerEvaluator<EvalT>(ev);
00260 }
00261 if (fieldManagerChoice == Albany::BUILD_RESID_FM) {
00262 PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter Stokes", dl->dummy);
00263 fm0.requireField<EvalT>(res_tag);
00264 }
00265 else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00266 fm0.template registerEvaluator<EvalT>
00267 (evalUtils.constructGatherSurfaceVelocityEvaluator());
00268 fm0.template registerEvaluator<EvalT>
00269 (evalUtils.constructGatherVelocityRMSEvaluator());
00270
00271 Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00272 return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
00273 }
00274
00275 return Teuchos::null;
00276 }
00277 #endif // FELIX_STOKESFOPROBLEM_HPP
