00001
00002
00003
00004
00005
00006 #include "Teuchos_TestForException.hpp"
00007 #include "Teuchos_CommHelpers.hpp"
00008
00009 namespace PHAL {
00010
00011
00012 void split(const std::string &s, char delim, std::vector<std::string> &elems) {
00013 std::stringstream ss(s);
00014 std::string item;
00015 while(std::getline(ss, item, delim)) {
00016 elems.push_back(item);
00017 }
00018 }
00019
00020 }
00021
00022 template<typename EvalT, typename Traits>
00023 PHAL::ResponseFieldIntegral<EvalT, Traits>::
00024 ResponseFieldIntegral(Teuchos::ParameterList& p,
00025 const Teuchos::RCP<Albany::Layouts>& dl) :
00026 coordVec("Coord Vec", dl->qp_gradient),
00027 weights("Weights", dl->qp_scalar)
00028 {
00029
00030 Teuchos::ParameterList* plist =
00031 p.get<Teuchos::ParameterList*>("Parameter List");
00032 Teuchos::RCP<const Teuchos::ParameterList> reflist =
00033 this->getValidResponseParameters();
00034 plist->validateParameters(*reflist,0);
00035
00036
00037 std::string field_name = plist->get<std::string>("Field Name");
00038 std::string fieldType = plist->get<std::string>("Field Type", "Scalar");
00039 if (plist->isType< Teuchos::Array<int> >("Field Components"))
00040 field_components = plist->get< Teuchos::Array<int> >("Field Components");
00041 Teuchos::RCP<PHX::DataLayout> field_layout;
00042 Teuchos::RCP<PHX::DataLayout> local_response_layout;
00043 Teuchos::RCP<PHX::DataLayout> global_response_layout;
00044 if (fieldType == "Scalar") {
00045 field_layout = dl->qp_scalar;
00046 local_response_layout = dl->cell_scalar;
00047 global_response_layout = dl->workset_scalar;
00048 }
00049 else if (fieldType == "Vector") {
00050 field_layout = dl->qp_vector;
00051 if (field_components.size() == 0) {
00052 local_response_layout = dl->cell_vector;
00053 global_response_layout = dl->workset_vector;
00054 }
00055 else {
00056 int worksetSize = dl->cell_scalar->dimension(0);
00057 local_response_layout =
00058 Teuchos::rcp(new PHX::MDALayout<Cell,Dim>(worksetSize,
00059 field_components.size()));
00060 global_response_layout =
00061 Teuchos::rcp(new PHX::MDALayout<Dim>(field_components.size()));
00062 }
00063 }
00064 else if (fieldType == "Tensor") {
00065 field_layout = dl->qp_tensor;
00066 local_response_layout = dl->cell_tensor;
00067 global_response_layout = dl->workset_tensor;
00068 }
00069 else {
00070 TEUCHOS_TEST_FOR_EXCEPTION(
00071 true,
00072 Teuchos::Exceptions::InvalidParameter,
00073 "Invalid field type " << fieldType << ". Support values are " <<
00074 "Scalar, Vector, and Tensor." << std::endl);
00075 }
00076 field = PHX::MDField<ScalarT>(field_name, field_layout);
00077 field_layout->dimensions(field_dims);
00078 field_rank = field_layout->rank();
00079 if (field_components.size() == 0) {
00080 int num_components = field_dims[field_rank-1];
00081 field_components.resize(num_components);
00082 for (int i=0; i<num_components; i++)
00083 field_components[i] = i;
00084 }
00085
00086
00087 std::vector<PHX::DataLayout::size_type> coord_dims;
00088 dl->qp_vector->dimensions(coord_dims);
00089 numQPs = coord_dims[1];
00090 numDims = coord_dims[2];
00091
00092
00093 std::string ebNameStr = plist->get<std::string>("Element Block Name","");
00094 if(ebNameStr.length() > 0) split(ebNameStr,',',ebNames);
00095
00096 limitX = limitY = limitZ = false;
00097 if( plist->isParameter("x min") && plist->isParameter("x max") ) {
00098 limitX = true; TEUCHOS_TEST_FOR_EXCEPT(numDims <= 0);
00099 xmin = plist->get<double>("x min");
00100 xmax = plist->get<double>("x max");
00101 }
00102 if( plist->isParameter("y min") && plist->isParameter("y max") ) {
00103 limitY = true; TEUCHOS_TEST_FOR_EXCEPT(numDims <= 1);
00104 ymin = plist->get<double>("y min");
00105 ymax = plist->get<double>("y max");
00106 }
00107 if( plist->isParameter("z min") && plist->isParameter("z max") ) {
00108 limitZ = true; TEUCHOS_TEST_FOR_EXCEPT(numDims <= 2);
00109 zmin = plist->get<double>("z min");
00110 zmax = plist->get<double>("z max");
00111 }
00112
00113
00114 double X0 = plist->get<double>("Length Scaling", 1.0);
00115 if (numDims == 1)
00116 scaling = X0;
00117 else if (numDims == 2)
00118 scaling = X0*X0;
00119 else if (numDims == 3)
00120 scaling = X0*X0*X0;
00121 else
00122 TEUCHOS_TEST_FOR_EXCEPTION(
00123 true, Teuchos::Exceptions::InvalidParameter,
00124 std::endl << "Error! Invalid number of dimensions: " << numDims <<
00125 std::endl);
00126
00127
00128 this->addDependentField(field);
00129 this->addDependentField(coordVec);
00130 this->addDependentField(weights);
00131 this->setName(field_name+" Response Field Integral"+PHX::TypeString<EvalT>::value);
00132
00133
00134 p.set("Stand-alone Evaluator", false);
00135 std::string local_response_name =
00136 field_name + " Local Response Field Integral";
00137 std::string global_response_name =
00138 field_name + " Global Response Field Integral";
00139 PHX::Tag<ScalarT> local_response_tag(local_response_name,
00140 local_response_layout);
00141 PHX::Tag<ScalarT> global_response_tag(global_response_name,
00142 global_response_layout);
00143 p.set("Local Response Field Tag", local_response_tag);
00144 p.set("Global Response Field Tag", global_response_tag);
00145 PHAL::SeparableScatterScalarResponse<EvalT,Traits>::setup(p,dl);
00146 }
00147
00148
00149 template<typename EvalT, typename Traits>
00150 void PHAL::ResponseFieldIntegral<EvalT, Traits>::
00151 postRegistrationSetup(typename Traits::SetupData d,
00152 PHX::FieldManager<Traits>& fm)
00153 {
00154 this->utils.setFieldData(field,fm);
00155 this->utils.setFieldData(coordVec,fm);
00156 this->utils.setFieldData(weights,fm);
00157 PHAL::SeparableScatterScalarResponse<EvalT,Traits>::postRegistrationSetup(d,fm);
00158 }
00159
00160
00161 template<typename EvalT, typename Traits>
00162 void PHAL::ResponseFieldIntegral<EvalT, Traits>::
00163 preEvaluate(typename Traits::PreEvalData workset)
00164 {
00165 for (typename PHX::MDField<ScalarT>::size_type i=0;
00166 i<this->global_response.size(); i++)
00167 this->global_response[i] = 0.0;
00168
00169
00170 PHAL::SeparableScatterScalarResponse<EvalT,Traits>::preEvaluate(workset);
00171 }
00172
00173
00174 template<typename EvalT, typename Traits>
00175 void PHAL::ResponseFieldIntegral<EvalT, Traits>::
00176 evaluateFields(typename Traits::EvalData workset)
00177 {
00178
00179 for (typename PHX::MDField<ScalarT>::size_type i=0;
00180 i<this->local_response.size(); i++)
00181 this->local_response[i] = 0.0;
00182
00183 if( ebNames.size() == 0 ||
00184 std::find(ebNames.begin(), ebNames.end(), workset.EBName) != ebNames.end() ) {
00185
00186 ScalarT s;
00187 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00188
00189 bool cellInBox = false;
00190 for (std::size_t qp=0; qp < numQPs; ++qp) {
00191 if( (!limitX || (coordVec(cell,qp,0) >= xmin && coordVec(cell,qp,0) <= xmax)) &&
00192 (!limitY || (coordVec(cell,qp,1) >= ymin && coordVec(cell,qp,1) <= ymax)) &&
00193 (!limitZ || (coordVec(cell,qp,2) >= zmin && coordVec(cell,qp,2) <= zmax)) ) {
00194 cellInBox = true; break; }
00195 }
00196 if( !cellInBox ) continue;
00197
00198 for (std::size_t qp=0; qp < numQPs; ++qp) {
00199 if (field_rank == 2) {
00200 s = field(cell,qp) * weights(cell,qp) * scaling;
00201 this->local_response(cell) += s;
00202 this->global_response(0) += s;
00203 }
00204 else if (field_rank == 3) {
00205 for (std::size_t dim=0; dim < field_components.size(); ++dim) {
00206 s = field(cell,qp,field_components[dim]) * weights(cell,qp) * scaling;
00207 this->local_response(cell,dim) += s;
00208 this->global_response(dim) += s;
00209 }
00210 }
00211 else if (field_rank == 4) {
00212 for (std::size_t dim1=0; dim1 < field_dims[2]; ++dim1) {
00213 for (std::size_t dim2=0; dim2 < field_dims[3]; ++dim2) {
00214 s = field(cell,qp,dim1,dim2) * weights(cell,qp) * scaling;
00215 this->local_response(cell,dim1,dim2) += s;
00216 this->global_response(dim1,dim2) += s;
00217 }
00218 }
00219 }
00220 }
00221 }
00222 }
00223
00224
00225 PHAL::SeparableScatterScalarResponse<EvalT,Traits>::evaluateFields(workset);
00226 }
00227
00228
00229 template<typename EvalT, typename Traits>
00230 void PHAL::ResponseFieldIntegral<EvalT, Traits>::
00231 postEvaluate(typename Traits::PostEvalData workset)
00232 {
00233
00234 Teuchos::RCP< Teuchos::ValueTypeSerializer<int,ScalarT> > serializer =
00235 workset.serializerManager.template getValue<EvalT>();
00236 Teuchos::reduceAll(
00237 *workset.comm, *serializer, Teuchos::REDUCE_SUM,
00238 this->global_response.size(), &this->global_response[0],
00239 &this->global_response[0]);
00240
00241
00242 PHAL::SeparableScatterScalarResponse<EvalT,Traits>::postEvaluate(workset);
00243 }
00244
00245
00246 template<typename EvalT,typename Traits>
00247 Teuchos::RCP<const Teuchos::ParameterList>
00248 PHAL::ResponseFieldIntegral<EvalT,Traits>::
00249 getValidResponseParameters() const
00250 {
00251 Teuchos::RCP<Teuchos::ParameterList> validPL =
00252 rcp(new Teuchos::ParameterList("Valid ResponseFieldIntegral Params"));
00253 Teuchos::RCP<const Teuchos::ParameterList> baseValidPL =
00254 PHAL::SeparableScatterScalarResponse<EvalT,Traits>::getValidResponseParameters();
00255 validPL->setParameters(*baseValidPL);
00256
00257 validPL->set<std::string>("Name", "", "Name of response function");
00258 validPL->set<int>("Phalanx Graph Visualization Detail", 0, "Make dot file to visualize phalanx graph");
00259 validPL->set<std::string>("Field Type", "", "Type of field (scalar, vector, ...)");
00260 validPL->set<std::string>(
00261 "Element Block Name", "",
00262 "Name of the element block to use as the integration domain");
00263 validPL->set<std::string>("Field Name", "", "Field to integrate");
00264 validPL->set<bool>("Positive Return Only",false);
00265
00266 validPL->set<double>("Length Scaling", 1.0, "Length Scaling");
00267 validPL->set<double>("x min", 0.0, "Integration domain minimum x coordinate");
00268 validPL->set<double>("x max", 0.0, "Integration domain maximum x coordinate");
00269 validPL->set<double>("y min", 0.0, "Integration domain minimum y coordinate");
00270 validPL->set<double>("y max", 0.0, "Integration domain maximum y coordinate");
00271 validPL->set<double>("z min", 0.0, "Integration domain minimum z coordinate");
00272 validPL->set<double>("z max", 0.0, "Integration domain maximum z coordinate");
00273
00274 validPL->set< Teuchos::Array<int> >("Field Components", Teuchos::Array<int>(),
00275 "Field components to scatter");
00276
00277 return validPL;
00278 }
00279
00280
00281
