• Main Page
  • Related Pages
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

Adapt_ElementSizeField_Def.hpp

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 
00007 #include <fstream>
00008 #include "Teuchos_TestForException.hpp"
00009 #include "Adapt_NodalDataBlock.hpp"
00010 
00011 template<typename T>
00012 T Sqr(T num)
00013 {
00014     return num * num;
00015 }
00016 
00017 template<typename EvalT, typename Traits>
00018 Adapt::ElementSizeFieldBase<EvalT, Traits>::
00019 ElementSizeFieldBase(Teuchos::ParameterList& p,
00020       const Teuchos::RCP<Albany::Layouts>& dl) :
00021   coordVec(p.get<std::string>("Coordinate Vector Name"), dl->qp_vector),
00022   coordVec_vertices(p.get<std::string>("Coordinate Vector Name"), dl->vertices_vector),
00023   qp_weights("Weights", dl->qp_scalar)
00024 {
00025 
00027   Teuchos::ParameterList* plist = 
00028     p.get<Teuchos::ParameterList*>("Parameter List");
00029   Teuchos::RCP<const Teuchos::ParameterList> reflist = 
00030     this->getValidSizeFieldParameters();
00031   plist->validateParameters(*reflist,0);
00032 
00033   // Isotropic --> element size scalar corresponding to the nominal element radius
00034   // Anisotropic --> element size vector (x, y, z) with the width, length, and height of the element
00035   // Weighted versions (upcoming) --> scale the above sizes with a scalar or vector field
00036 
00037   className = plist->get<std::string>("Size Field Name", "Element_Size_Field");
00038   outputToExodus = plist->get<bool>("Output to File", true);
00039   outputCellAverage = plist->get<bool>("Generate Cell Average", true);
00040   outputQPData = plist->get<bool>("Generate QP Values", false);
00041   outputNodeData = plist->get<bool>("Generate Nodal Values", false);
00042   isAnisotropic = plist->get<bool>("Anisotropic Size Field", false);
00043 
00045   Teuchos::RCP<PHX::DataLayout> scalar_dl = dl->qp_scalar;
00046   Teuchos::RCP<PHX::DataLayout> vector_dl = dl->qp_vector;
00047   Teuchos::RCP<PHX::DataLayout> cell_dl = dl->cell_scalar;
00048   Teuchos::RCP<PHX::DataLayout> vert_vector_dl = dl->vertices_vector;
00049   numQPs = vector_dl->dimension(1);
00050   numDims = vector_dl->dimension(2);
00051   numVertices = vert_vector_dl->dimension(2);
00052  
00053   this->addDependentField(qp_weights);
00054   this->addDependentField(coordVec);
00055   this->addDependentField(coordVec_vertices);
00056 
00058   this->pStateMgr = p.get< Albany::StateManager* >("State Manager Ptr");
00059 
00060   if( outputCellAverage ) {
00061     if(isAnisotropic) //An-isotropic
00062       this->pStateMgr->registerStateVariable(className + "_Cell", dl->cell_vector, dl->dummy, "all", "scalar", 
00063          0.0, false, outputToExodus);
00064     else
00065       this->pStateMgr->registerStateVariable(className + "_Cell", dl->cell_scalar, dl->dummy, "all", "scalar", 
00066          0.0, false, outputToExodus);
00067   }
00068 
00069   if( outputQPData ) {
00070 //    if(isAnisotropic) //An-isotropic
00071 //    Always anisotropic?
00072       this->pStateMgr->registerStateVariable(className + "_QP", dl->qp_vector, dl->dummy, "all", 
00073         "scalar", 0.0, false, outputToExodus);
00074 //    else
00075 //      this->pStateMgr->registerStateVariable(className + "_QP", dl->qp_scalar, dl->dummy, "all", 
00076 //        "scalar", 0.0, false, outputToExodus);
00077   }
00078 
00079   if( outputNodeData ) {
00080     // The weighted projected value
00081 
00082     // Note that all dl->node_node_* layouts are handled by the Adapt_NodalDataBlock class, inside
00083     // of the state manager, as they require interprocessor synchronization
00084 
00085     if(isAnisotropic){ //An-isotropic
00086       this->pStateMgr->registerStateVariable(className + "_Node", dl->node_node_vector, dl->dummy, "all", 
00087          "scalar", 0.0, false, outputToExodus);
00088 
00089 
00090     }
00091     else {
00092       this->pStateMgr->registerStateVariable(className + "_Node", dl->node_node_scalar, dl->dummy, "all", 
00093          "scalar", 0.0, false, outputToExodus);
00094 
00095     }
00096 
00097     // The value of the weights used in the projection
00098     // Initialize to zero - should give us nan's during the division step if something is wrong
00099     this->pStateMgr->registerStateVariable(className + "_NodeWgt", dl->node_node_scalar, dl->dummy, "all", 
00100          "scalar", 0.0, false, outputToExodus);
00101 
00102   }
00103 
00104   // Create field tag
00105   size_field_tag = 
00106     Teuchos::rcp(new PHX::Tag<ScalarT>(className, dl->dummy));
00107 
00108   this->addEvaluatedField(*size_field_tag);
00109 }
00110 
00111 // **********************************************************************
00112 template<typename EvalT, typename Traits>
00113 void Adapt::ElementSizeFieldBase<EvalT, Traits>::
00114 postRegistrationSetup(typename Traits::SetupData d,
00115                       PHX::FieldManager<Traits>& fm)
00116 {
00117 //  this->utils.setFieldData(field,fm);
00118   this->utils.setFieldData(qp_weights,fm);
00119   this->utils.setFieldData(coordVec,fm);
00120   this->utils.setFieldData(coordVec_vertices,fm);
00121 
00122 }
00123 
00124 // **********************************************************************
00125 // Specialization: Residual
00126 // **********************************************************************
00127 // **********************************************************************
00128 template<typename Traits>
00129 Adapt::
00130 ElementSizeField<PHAL::AlbanyTraits::Residual, Traits>::
00131 ElementSizeField(Teuchos::ParameterList& p, const Teuchos::RCP<Albany::Layouts>& dl) :
00132   ElementSizeFieldBase<PHAL::AlbanyTraits::Residual, Traits>(p, dl)
00133 {
00134 }
00135 
00136 template<typename Traits>
00137 void Adapt::ElementSizeField<PHAL::AlbanyTraits::Residual, Traits>::
00138 preEvaluate(typename Traits::PreEvalData workset)
00139 {
00140   // Note that we only need to initialize the vectors when dealing with node data, as we assume
00141   // the vectors are initialized to zero for Epetra_Export "ADD" operation
00142   // Zero data for accumulation here
00143   if( this->outputNodeData ) { 
00144     Teuchos::RCP<Adapt::NodalDataBlock> node_data = this->pStateMgr->getStateInfoStruct()->getNodalDataBlock();
00145     node_data->initializeVectors(0.0);
00146   }
00147 }
00148 
00149 template<typename Traits>
00150 void Adapt::ElementSizeField<PHAL::AlbanyTraits::Residual, Traits>::
00151 evaluateFields(typename Traits::EvalData workset)
00152 {
00153   double value;
00154 
00155   if( this->outputCellAverage ) { // nominal radius
00156 
00157     // Get shards Array (from STK) for this workset
00158     Albany::MDArray data = (*workset.stateArrayPtr)[this->className + "_Cell"];
00159     std::vector<int> dims;
00160     data.dimensions(dims);
00161     int size = dims.size();
00162 
00163 
00164     for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00165           this->getCellRadius(cell, value);
00166 //          data(cell, (std::size_t)0) = ADValue(value);
00167           data(cell, (std::size_t)0) = value;
00168     }
00169   }
00170 
00171   if( this->outputQPData ) { // x_\xi \cdot x_\xi, x_\eta \cdot x_\eta, x_\zeta \cdot x_\zeta
00172 
00173     // Get shards Array (from STK) for this workset
00174     Albany::MDArray data = (*workset.stateArrayPtr)[this->className + "_QP"];
00175     std::vector<int> dims;
00176     data.dimensions(dims);
00177     int size = dims.size();
00178 
00179 
00180     for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00181           this->getCellRadius(cell, value);
00182 //          data(cell, (std::size_t)0) = ADValue(value);
00183           data(cell, (std::size_t)0) = value;
00184     }
00185 /*
00186     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00187       for (std::size_t qp=0; qp < numQPs; ++qp) {      
00188         for (std::size_t i=0; i < numDims; ++i) { // loop over \xi, \eta, \zeta
00189           data(cell, qp, i) = 0.0;
00190           for (std::size_t j=0; j < numDims; ++j) {
00191             data(cell, qp, i) += coordVec(cell, qp, j) * wGradBF(cell, node, qp, j);
00192             for (std::size_t alpha=0; alpha < numDims; ++alpha) {  
00193               Gc(cell,qp,i,j) += jacobian_inv(cell,qp,alpha,i)*jacobian_inv(cell,qp,alpha,j); 
00194             }
00195           } 
00196         } 
00197       }
00198     }
00199 */
00200 
00201   }
00202 
00203   if( this->outputNodeData ) { // nominal radius, store as nodal data that will be scattered and summed
00204 
00205     // Get the node data block container
00206     Teuchos::RCP<Adapt::NodalDataBlock> node_data = this->pStateMgr->getStateInfoStruct()->getNodalDataBlock();
00207     Teuchos::RCP<Epetra_Vector> data = node_data->getLocalNodeVec();
00208     Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >  wsElNodeID = workset.wsElNodeID;
00209     Teuchos::RCP<const Epetra_BlockMap> local_node_map = node_data->getLocalMap();
00210 
00211     int l_nV = this->numVertices;
00212     int l_nD = this->numDims;
00213     int blocksize = node_data->getBlocksize();
00214 
00215     int  node_var_offset;
00216     int  node_var_ndofs;
00217     int  node_weight_offset;
00218     int  node_weight_ndofs;
00219     node_data->getNDofsAndOffset(this->className + "_Node", node_var_offset, node_var_ndofs);
00220     node_data->getNDofsAndOffset(this->className + "_NodeWgt", node_weight_offset, node_weight_ndofs);
00221 
00222     for (int cell = 0; cell < workset.numCells; ++cell) { // loop over all elements in workset
00223 
00224       std::vector<double> maxCoord(3,-1e10);
00225       std::vector<double> minCoord(3,+1e10);
00226 
00227       // Get element width in x, y, z
00228       for (int v=0; v < l_nV; ++v) { // loop over all the "corners" of each element
00229         for (int k=0; k < l_nD; ++k) { // loop over each dimension of the problem
00230           if(maxCoord[k] < this->coordVec_vertices(cell,v,k)) maxCoord[k] = this->coordVec_vertices(cell,v,k);
00231           if(minCoord[k] > this->coordVec_vertices(cell,v,k)) minCoord[k] = this->coordVec_vertices(cell,v,k);
00232         }
00233       }
00234 
00235       if(this->isAnisotropic) //An-isotropic
00236         // Note: code assumes blocksize of blockmap is numDims + 1 - the last entry accumulates the weight
00237         for (int node = 0; node < l_nV; ++node) { // loop over all the "corners" of each element
00238 
00239           int global_node = wsElNodeID[cell][node]; // get the global id of this node
00240 
00241           int local_node = local_node_map->LID(global_node); // skip the node if it is not owned by me
00242           if(local_node < 0) continue;
00243 
00244           // accumulate 1/2 of the element width in each dimension - into each element corner
00245           for (int k=0; k < node_var_ndofs; ++k) 
00246 //            data[global_node][k] += ADValue(maxCoord[k] - minCoord[k]) / 2.0;
00247             (*data)[local_node * blocksize + node_var_offset + k] += (maxCoord[k] - minCoord[k]) / 2.0;
00248 
00249           // save the weight (denominator)
00250           (*data)[local_node * blocksize + node_weight_offset] += 1.0;
00251 
00252       } // end anisotropic size field
00253 
00254       else // isotropic size field
00255         // Note: code assumes blocksize of blockmap is 1 + 1 = 2 - the last entry accumulates the weight
00256         for (int node = 0; node < l_nV; ++node) { // loop over all the "corners" of each element
00257 
00258           int global_node = wsElNodeID[cell][node]; // get the global id of this node
00259 
00260           int local_node = local_node_map->LID(global_node); // skip the node if it is not owned by me
00261           if(local_node < 0) continue;
00262 
00263           // save element radius, just a scalar
00264           for (int k=0; k < l_nD; ++k) {
00265 //            data[global_node][k] += ADValue(maxCoord[k] - minCoord[k]) / 2.0;
00266             (*data)[local_node * blocksize + node_var_offset] += (maxCoord[k] - minCoord[k]) / 2.0;
00267             // save the weight (denominator)
00268             (*data)[local_node * blocksize + node_weight_offset] += 1.0;
00269 
00270           }
00271 
00272 
00273           // the above calculates the average of the element width, depth, and height when
00274           // divided by the accumulated weights
00275 
00276       } // end isotropic size field
00277     } // end cell loop
00278   } // end node data if
00279 
00280 }
00281 
00282 template<typename Traits>
00283 void Adapt::ElementSizeField<PHAL::AlbanyTraits::Residual, Traits>::
00284 postEvaluate(typename Traits::PostEvalData workset)
00285 {
00286 
00287   if( this->outputNodeData ) {
00288 
00289     // Note: we are in postEvaluate so all PEs call this
00290 
00291     // Get the node data block container
00292     Teuchos::RCP<Adapt::NodalDataBlock> node_data = this->pStateMgr->getStateInfoStruct()->getNodalDataBlock();
00293     Teuchos::RCP<Epetra_Vector> data = node_data->getOverlapNodeVec();
00294     Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >  wsElNodeID = workset.wsElNodeID;
00295     Teuchos::RCP<const Epetra_BlockMap> overlap_node_map = node_data->getOverlapMap();
00296 
00297     int  node_var_offset;
00298     int  node_var_ndofs;
00299     int  node_weight_offset;
00300     int  node_weight_ndofs;
00301     node_data->getNDofsAndOffset(this->className + "_Node", node_var_offset, node_var_ndofs);
00302     node_data->getNDofsAndOffset(this->className + "_NodeWgt", node_weight_offset, node_weight_ndofs);
00303 
00304     // Build the exporter
00305     node_data->initializeExport();
00306 
00307     // do the export
00308     node_data->exportAddNodalDataBlock();
00309 
00310     int numNodes = overlap_node_map->NumMyElements();
00311     int blocksize = node_data->getBlocksize();
00312 
00313     // if isotropic, blocksize == 2 , if anisotropic blocksize == nDOF at node + 1
00314     // ndim if vector, ndim * ndim if tensor
00315 
00316     // all PEs divide the accumulated value(s) by the weights
00317 
00318     for (int overlap_node=0; overlap_node < numNodes; ++overlap_node)
00319 
00320       for (int k=0; k < node_var_ndofs; ++k) 
00321             (*data)[overlap_node * blocksize + node_var_offset + k] /=
00322                 (*data)[overlap_node * blocksize + node_weight_offset];
00323 
00324 
00325     // Export the data from the local to overlapped decomposition
00326     // Divide the overlap field through by the weights
00327     // Store the overlapped vector data back in stk in the field "field_name"
00328 
00329     node_data->saveNodalDataState();
00330 
00331   }
00332 
00333 }
00334 
00335 // **********************************************************************
00336 
00337 template<typename EvalT, typename Traits>
00338 void Adapt::ElementSizeFieldBase<EvalT, Traits>::
00339 getCellRadius(const std::size_t cell, typename EvalT::MeshScalarT& cellRadius) const
00340 {
00341   std::vector<MeshScalarT> maxCoord(3,-1e10);
00342   std::vector<MeshScalarT> minCoord(3,+1e10);
00343 
00344   for (int v=0; v < numVertices; ++v) {
00345     for (int k=0; k < numDims; ++k) {
00346       if(maxCoord[k] < coordVec_vertices(cell,v,k)) maxCoord[k] = coordVec_vertices(cell,v,k);
00347       if(minCoord[k] > coordVec_vertices(cell,v,k)) minCoord[k] = coordVec_vertices(cell,v,k);
00348     }
00349   }
00350 
00351   cellRadius = 0.0;
00352   for (int k=0; k < numDims; ++k) 
00353     cellRadius += (maxCoord[k] - minCoord[k]) *  (maxCoord[k] - minCoord[k]);
00354 
00355   cellRadius = std::sqrt(cellRadius) / 2.0;
00356 
00357 }
00358 
00359 
00360 // **********************************************************************
00361 template<typename EvalT, typename Traits>
00362 Teuchos::RCP<const Teuchos::ParameterList>
00363 Adapt::ElementSizeFieldBase<EvalT,Traits>::getValidSizeFieldParameters() const
00364 {
00365   Teuchos::RCP<Teuchos::ParameterList> validPL =
00366       rcp(new Teuchos::ParameterList("Valid ElementSizeField Params"));;
00367 
00368   validPL->set<std::string>("Name", "", "Name of size field Evaluator");
00369   validPL->set<std::string>("Size Field Name", "", "Size field prefix");
00370 
00371   validPL->set<bool>("Output to File", true, "Whether size field info should be output to a file");
00372   validPL->set<bool>("Generate Cell Average", true, "Whether cell average field should be generated");
00373   validPL->set<bool>("Generate QP Values", true, "Whether values at the quadpoints should be generated");
00374   validPL->set<bool>("Generate Nodal Values", true, "Whether values at the nodes should be generated");
00375   validPL->set<bool>("Anisotropic Size Field", true, "Is this size field calculation anisotropic?");
00376 
00377   return validPL;
00378 }
00379 
00380 // **********************************************************************
00381 

Generated on Wed Mar 26 2014 18:36:35 for Albany: a Trilinos-based PDE code by  doxygen 1.7.1