00001
00002
00003
00004
00005
00006
00007 #ifndef ALBANY_HELMHOLTZN2DPROBLEM_HPP
00008 #define ALBANY_HELMHOLTZN2DPROBLEM_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 namespace Albany {
00020
00025 class Helmholtz2DProblem : public Albany::AbstractProblem {
00026 public:
00027
00029 Helmholtz2DProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00030 const Teuchos::RCP<ParamLib>& paramLib);
00031
00033 virtual ~Helmholtz2DProblem();
00034
00036 virtual int spatialDimension() const { return 2; }
00037
00039 virtual void buildProblem(
00040 Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> > meshSpecs,
00041 StateManager& stateMgr);
00042
00043
00044 virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00045 buildEvaluators(
00046 PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00047 const Albany::MeshSpecsStruct& meshSpecs,
00048 Albany::StateManager& stateMgr,
00049 Albany::FieldManagerChoice fmchoice,
00050 const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00051
00053 Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00054
00055 private:
00056
00058 Helmholtz2DProblem(const Helmholtz2DProblem&);
00059
00061 Helmholtz2DProblem& operator=(const Helmholtz2DProblem&);
00062
00063 public:
00064
00066 template <typename EvalT> Teuchos::RCP<const PHX::FieldTag>
00067 constructEvaluators(
00068 PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00069 const Albany::MeshSpecsStruct& meshSpecs,
00070 Albany::StateManager& stateMgr,
00071 Albany::FieldManagerChoice fmchoice,
00072 const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00073
00074 void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00075
00076 protected:
00077
00079 double ksqr;
00080 bool haveSource;
00081 };
00082
00083 }
00084
00085 #include "Albany_Utils.hpp"
00086 #include "Albany_ProblemUtils.hpp"
00087 #include "Albany_EvaluatorUtils.hpp"
00088 #include "Albany_ResponseUtilities.hpp"
00089
00090 #include "PHAL_Source.hpp"
00091 #include "PHAL_HelmholtzResid.hpp"
00092
00093
00094 template <typename EvalT>
00095 Teuchos::RCP<const PHX::FieldTag>
00096 Albany::Helmholtz2DProblem::constructEvaluators(
00097 PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00098 const Albany::MeshSpecsStruct& meshSpecs,
00099 Albany::StateManager& stateMgr,
00100 Albany::FieldManagerChoice fieldManagerChoice,
00101 const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00102 {
00103 using Teuchos::RCP;
00104 using Teuchos::rcp;
00105 using Teuchos::ParameterList;
00106 using PHX::DataLayout;
00107 using PHX::MDALayout;
00108 using std::vector;
00109 using std::string;
00110 using PHAL::AlbanyTraits;
00111
00112 RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology(&meshSpecs.ctd));
00113 RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00114 intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00115
00116 const int numNodes = intrepidBasis->getCardinality();
00117 const int worksetSize = meshSpecs.worksetSize;
00118
00119 Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00120
00121 RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00122
00123 const int numDim = cubature->getDimension();
00124 const int numQPts = cubature->getNumPoints();
00125 const int numVertices = cellType->getNodeCount();
00126
00127 *out << "Field Dimensions: Workset=" << worksetSize
00128 << ", Vertices= " << numVertices
00129 << ", Nodes= " << numNodes
00130 << ", QuadPts= " << numQPts
00131 << ", Dim= " << numDim << std::endl;
00132
00133 RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00134 Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00135 bool supportsTransient=false;
00136
00137
00138 Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00139
00140
00141
00142 Teuchos::ArrayRCP<string> dof_names(neq);
00143 dof_names[0] = "U";
00144 dof_names[1] = "V";
00145
00146 Teuchos::ArrayRCP<string> dof_names_dot(neq);
00147 if (supportsTransient) {
00148 for (int i=0; i<neq; i++) dof_names_dot[i] = dof_names[i]+"_dot";
00149 }
00150
00151 Teuchos::ArrayRCP<string> resid_names(neq);
00152 for (int i=0; i<neq; i++) resid_names[i] = dof_names[i]+" Residual";
00153
00154 if (supportsTransient) fm0.template registerEvaluator<EvalT>
00155 (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot));
00156 else fm0.template registerEvaluator<EvalT>
00157 (evalUtils.constructGatherSolutionEvaluator_noTransient(false, dof_names));
00158
00159 fm0.template registerEvaluator<EvalT>
00160 (evalUtils.constructScatterResidualEvaluator(false, resid_names));
00161
00162 fm0.template registerEvaluator<EvalT>
00163 (evalUtils.constructGatherCoordinateVectorEvaluator());
00164
00165 fm0.template registerEvaluator<EvalT>
00166 (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00167
00168 fm0.template registerEvaluator<EvalT>
00169 (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00170
00171 for (int i=0; i<neq; i++) {
00172 fm0.template registerEvaluator<EvalT>
00173 (evalUtils.constructDOFInterpolationEvaluator(dof_names[i], i));
00174
00175 if (supportsTransient)
00176 fm0.template registerEvaluator<EvalT>
00177 (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[i], i));
00178
00179 fm0.template registerEvaluator<EvalT>
00180 (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i], i));
00181 }
00182
00183 if (haveSource) {
00184 RCP<ParameterList> p = rcp(new ParameterList("Helmholtz Source U"));
00185
00186 p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00187 p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00188
00189 p->set<string>("Pressure Source Name", "U GaussMonotone");
00190 p->set<string>("QP Coordinate Vector Name", "Coord Vec");
00191
00192
00193 p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00194 Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00195 p->set<Teuchos::ParameterList*>("Parameter List", ¶mList);
00196
00197 ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00198 fm0.template registerEvaluator<EvalT>(ev);
00199 }
00200
00201 if (haveSource) {
00202 RCP<ParameterList> p = rcp(new ParameterList("Helmholtz Source V"));
00203
00204 p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00205 p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00206
00207 p->set<string>("Pressure Source Name", "V GaussMonotone");
00208 p->set<string>("QP Coordinate Vector Name", "Coord Vec");
00209
00210
00211 p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00212 Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00213 p->set<Teuchos::ParameterList*>("Parameter List", ¶mList);
00214
00215 ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00216 fm0.template registerEvaluator<EvalT>(ev);
00217 }
00218
00219 {
00220 RCP<ParameterList> p = rcp(new ParameterList("Helmholtz Resid"));
00221
00222
00223 p->set<string>("Weighted BF Name", "wBF");
00224 p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00225
00226 p->set<string>("U Variable Name", "U");
00227 p->set<string>("V Variable Name", "V");
00228 p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00229
00230 p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00231 p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00232
00233 p->set<string>("U Gradient Variable Name", "U Gradient");
00234 p->set<string>("V Gradient Variable Name", "V Gradient");
00235 p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00236
00237 p->set<bool>("Have Source", haveSource);
00238 p->set<string>("U Pressure Source Name", "U GaussMonotone");
00239 p->set<string>("V Pressure Source Name", "V GaussMonotone");
00240
00241 p->set<double>("Ksqr", ksqr);
00242 p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00243
00244
00245 p->set<string>("U Residual Name", "U Residual");
00246 p->set<string>("V Residual Name", "V Residual");
00247 p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00248
00249 ev = rcp(new PHAL::HelmholtzResid<EvalT,AlbanyTraits>(*p));
00250 fm0.template registerEvaluator<EvalT>(ev);
00251 }
00252
00253
00254 if (fieldManagerChoice == Albany::BUILD_RESID_FM) {
00255 PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00256 fm0.requireField<EvalT>(res_tag);
00257 return res_tag.clone();
00258 }
00259
00260 else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM){
00261 Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00262 return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
00263 }
00264
00265 return Teuchos::null;
00266 }
00267 #endif // ALBANY_HEATNONLINEARSOURCE2DPROBLEM_HPP
