00001
00002
00003
00004
00005
00006
00007 #include <fstream>
00008 #include "Teuchos_TestForException.hpp"
00009
00010
00011 template<typename EvalT, typename Traits>
00012 QCAD::ResponseSaveField<EvalT, Traits>::
00013 ResponseSaveField(Teuchos::ParameterList& p,
00014 const Teuchos::RCP<Albany::Layouts>& dl) :
00015 weights("Weights", dl->qp_scalar)
00016 {
00018 Teuchos::ParameterList* plist =
00019 p.get<Teuchos::ParameterList*>("Parameter List");
00020 Teuchos::RCP<const Teuchos::ParameterList> reflist =
00021 this->getValidResponseParameters();
00022 plist->validateParameters(*reflist,0);
00023
00025 if(plist->isParameter("Vector Field Name")) {
00026 fieldName = plist->get<std::string>("Vector Field Name");
00027 isVectorField = true;
00028 }
00029 else {
00030 fieldName = plist->get<std::string>("Field Name");
00031 isVectorField = false;
00032 }
00033 stateName = plist->get<std::string>("State Name", fieldName);
00034 outputToExodus = plist->get<bool>("Output to Exodus", true);
00035 outputCellAverage = plist->get<bool>("Output Cell Average", true);
00036 memoryHolderOnly = plist->get<bool>("Memory Placeholder Only", false);
00037 vectorOp = plist->get<std::string>("Vector Operation", "magnitude");
00038
00040 Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
00041 Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;
00042 Teuchos::RCP<PHX::DataLayout> cell_dl = dl->cell_scalar;
00043 numQPs = vector_dl->dimension(1);
00044 numDims = vector_dl->dimension(2);
00045
00047 Teuchos::RCP<PHX::DataLayout>& field_dl = isVectorField ? vector_dl : scalar_dl;
00048 PHX::MDField<ScalarT> f(fieldName, field_dl); field = f;
00049 this->addDependentField(field);
00050 this->addDependentField(weights);
00051
00053 Albany::StateManager* pStateMgr = p.get< Albany::StateManager* >("State Manager Ptr");
00054 if( outputCellAverage ) {
00055 pStateMgr->registerStateVariable(stateName, cell_dl, "ALL", "scalar", 0.0, false, outputToExodus);
00056 }
00057 else {
00058 pStateMgr->registerStateVariable(stateName, scalar_dl, "ALL", "scalar", 0.0, false, outputToExodus);
00059 }
00060
00061
00062 response_field_tag =
00063 Teuchos::rcp(new PHX::Tag<ScalarT>(fieldName + " Save Field Response",
00064 dl->dummy));
00065 this->addEvaluatedField(*response_field_tag);
00066 }
00067
00068
00069 template<typename EvalT, typename Traits>
00070 void QCAD::ResponseSaveField<EvalT, Traits>::
00071 postRegistrationSetup(typename Traits::SetupData d,
00072 PHX::FieldManager<Traits>& fm)
00073 {
00074 this->utils.setFieldData(field,fm);
00075 this->utils.setFieldData(weights,fm);
00076 }
00077
00078
00079 template<typename EvalT, typename Traits>
00080 void QCAD::ResponseSaveField<EvalT, Traits>::
00081 evaluateFields(typename Traits::EvalData workset)
00082 {
00083 using Albany::ADValue;
00084
00085 const std::size_t iX=0;
00086 const std::size_t iY=1;
00087 const std::size_t iZ=2;
00088
00089
00090
00091 if(memoryHolderOnly) return;
00092
00093
00094
00095 Albany::MDArray sta = (*workset.stateArrayPtr)[stateName];
00096 std::vector<int> dims;
00097 sta.dimensions(dims);
00098 int size = dims.size();
00099
00100 if(!isVectorField) {
00101 switch (size) {
00102 case 2:
00103 for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00104 if( outputCellAverage ) {
00105 double integral = 0, vol = 0;
00106 for (std::size_t qp = 0; qp < numQPs; ++qp) {
00107 integral += ADValue(field(cell,qp)) * ADValue(weights(cell,qp));
00108 vol += ADValue(weights(cell, qp));
00109 }
00110 sta(cell,(std::size_t)0) = integral / vol;
00111 }
00112 else {
00113 for (std::size_t qp = 0; qp < numQPs; ++qp)
00114 sta(cell, qp) = ADValue(field(cell,qp));
00115 }
00116 }
00117 break;
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132 default:
00133 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00134 "Unexpected dimensions in SaveField response Evaluator: " << size);
00135 }
00136 }
00137 else {
00138 ScalarT t;
00139 switch (size) {
00140 case 2:
00141 for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00142
00143 ScalarT stateValue = 0.0;
00144 double vol = 0.0;
00145 if( outputCellAverage ) sta(cell,(std::size_t)0) = 0.0;
00146
00147 for (std::size_t qp = 0; qp < numQPs; ++qp) {
00148 t = 0.0;
00149
00150 if(vectorOp == "magnitude") {
00151 for (std::size_t i = 0; i < numDims; ++i)
00152 t += field(cell,qp,i)*field(cell,qp,i);
00153 stateValue = sqrt(t);
00154 }
00155 else if(vectorOp == "xyMagnitude") {
00156 if(numDims > iX) t += field(cell,qp,iX)*field(cell,qp,iX);
00157 if(numDims > iY) t += field(cell,qp,iY)*field(cell,qp,iY);
00158 stateValue = sqrt(t);
00159 }
00160 else if(vectorOp == "xzMagnitude") {
00161 if(numDims > iX) t += field(cell,qp,iX)*field(cell,qp,iX);
00162 if(numDims > iZ) t += field(cell,qp,iZ)*field(cell,qp,iZ);
00163 stateValue = sqrt(t);
00164 }
00165 else if(vectorOp == "yzMagnitude") {
00166 if(numDims > iY) t += field(cell,qp,iY)*field(cell,qp,iY);
00167 if(numDims > iZ) t += field(cell,qp,iZ)*field(cell,qp,iZ);
00168 stateValue = sqrt(t);
00169 }
00170
00171 else if(vectorOp == "magnitude2") {
00172 for (std::size_t i = 0; i < numDims; ++i)
00173 t += field(cell,qp,i)*field(cell,qp,i);
00174 stateValue = t;
00175 }
00176 else if(vectorOp == "xyMagnitude2") {
00177 if(numDims > iX) t += field(cell,qp,iX)*field(cell,qp,iX);
00178 if(numDims > iY) t += field(cell,qp,iY)*field(cell,qp,iY);
00179 stateValue = t;
00180 }
00181 else if(vectorOp == "xzMagnitude2") {
00182 if(numDims > iX) t += field(cell,qp,iX)*field(cell,qp,iX);
00183 if(numDims > iZ) t += field(cell,qp,iZ)*field(cell,qp,iZ);
00184 stateValue = t;
00185 }
00186 else if(vectorOp == "yzMagnitude2") {
00187 if(numDims > iY) t += field(cell,qp,iY)*field(cell,qp,iY);
00188 if(numDims > iZ) t += field(cell,qp,iZ)*field(cell,qp,iZ);
00189 stateValue = t;
00190 }
00191
00192 else if(vectorOp == "xCoord") {
00193 if(numDims > iX) stateValue = field(cell,qp,iX);
00194 }
00195 else if(vectorOp == "yCoord") {
00196 if(numDims > iY) stateValue = field(cell,qp,iY);
00197 }
00198 else if(vectorOp == "zCoord") {
00199 if(numDims > iZ) stateValue = field(cell,qp,iZ);
00200 }
00201 else {
00202 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00203 "Unknown vector operation: " << vectorOp);
00204 }
00205
00206 if( outputCellAverage ) {
00207 sta(cell, (std::size_t)0) += ADValue(stateValue) * ADValue(weights(cell,qp));
00208 vol += ADValue(weights(cell,qp));
00209 }
00210 else sta(cell, qp) = ADValue(stateValue);
00211 }
00212
00213 if( outputCellAverage ) sta(cell,(std::size_t)0) /= vol;
00214 }
00215 break;
00216 default:
00217 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00218 "Unexpected dimensions in SaveField response Evaluator: " << size);
00219 }
00220 }
00221 }
00222
00223
00224 template<typename EvalT, typename Traits>
00225 Teuchos::RCP<const Teuchos::ParameterList>
00226 QCAD::ResponseSaveField<EvalT,Traits>::getValidResponseParameters() const
00227 {
00228 Teuchos::RCP<Teuchos::ParameterList> validPL =
00229 rcp(new Teuchos::ParameterList("Valid ResponseSaveField Params"));;
00230
00231 validPL->set<std::string>("Name", "", "Name of response function");
00232 validPL->set<int>("Phalanx Graph Visualization Detail", 0, "Make dot file to visualize phalanx graph");
00233 validPL->set<std::string>("Field Name", "", "Field to save");
00234 validPL->set<std::string>("Vector Field Name", "", "Vector field to save");
00235 validPL->set<std::string>("Vector Operation", "magnitude", "How to convert vector to scalar value, e.g., magnitude, xyMagnitude, xCoord");
00236 validPL->set<std::string>("State Name", "<Field Name>", "State name to save field as");
00237 validPL->set<bool>("Output to Exodus", true, "Whether state should be output in STK dump to exodus");
00238 validPL->set<bool>("Output Cell Average", true, "Whether cell average or all quadpoint data should be output to exodus");
00239 validPL->set<bool>("Memory Placeholder Only", false, "True if data should not actually be transferred to this state, i.e., the state is just used as a memory container and should not be overwritten when responses are computed");
00240 validPL->set<std::string>("Description", "", "Description of this response used by post processors");
00241
00242 return validPL;
00243 }
00244
00245
00246
