Albany: a Trilinos-based PDE code: Albany_STKDiscretization.hpp Source File

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 #ifndef ALBANY_STKDISCRETIZATION_HPP
00008 #define ALBANY_STKDISCRETIZATION_HPP
00009 
00010 #include <vector>
00011 #include <utility>
00012 
00013 #include "Teuchos_ParameterList.hpp"
00014 #include "Teuchos_VerboseObject.hpp"
00015 
00016 #include "Epetra_Comm.h"
00017 
00018 #include "Albany_AbstractDiscretization.hpp"
00019 #include "Albany_AbstractSTKMeshStruct.hpp"
00020 
00021 #include "Epetra_CrsMatrix.h"
00022 #include "Epetra_Vector.h"
00023 
00024 #include "Piro_NullSpaceUtils.hpp" // has defn of struct that holds null space info for ML
00025 
00026 // Start of STK stuff
00027 #include <stk_util/parallel/Parallel.hpp>
00028 #include <stk_mesh/base/Types.hpp>
00029 #include <stk_mesh/base/MetaData.hpp>
00030 #include <stk_mesh/base/BulkData.hpp>
00031 #include <stk_mesh/base/Field.hpp>
00032 #include <stk_mesh/base/FieldTraits.hpp>
00033 #ifdef ALBANY_SEACAS
00034   #include <stk_io/MeshReadWriteUtils.hpp>
00035 #endif
00036 
00037 
00038 namespace Albany {
00039 
00040   struct MeshGraph {
00041 
00042        std::vector<std::size_t> start;
00043        std::vector<std::size_t> adj;
00044 
00045   };
00046 
00047   class STKDiscretization : public Albany::AbstractDiscretization {
00048   public:
00049 
00051     STKDiscretization(
00052        Teuchos::RCP<Albany::AbstractSTKMeshStruct> stkMeshStruct,
00053        const Teuchos::RCP<const Epetra_Comm>& comm,
00054        const Teuchos::RCP<Piro::MLRigidBodyModes>& rigidBodyModes = Teuchos::null);
00055 
00056 
00058     ~STKDiscretization();
00059 
00061     Teuchos::RCP<const Epetra_Map> getMap() const;
00062 
00064     Teuchos::RCP<const Epetra_Map> getOverlapMap() const;
00065 
00067     Teuchos::RCP<const Epetra_CrsGraph> getJacobianGraph() const;
00068 
00070     Teuchos::RCP<const Epetra_CrsGraph> getOverlapJacobianGraph() const;
00071 
00073     Teuchos::RCP<const Epetra_Map> getNodeMap() const;
00074 
00076     const NodeSetList& getNodeSets() const { return nodeSets; };
00077     const NodeSetCoordList& getNodeSetCoords() const { return nodeSetCoords; };
00078 
00080     const SideSetList& getSideSets(const int workset) const { return sideSets[workset]; };
00081 
00083     WsLIDList& getElemGIDws() { return elemGIDws; };
00084 
00086     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type& getWsElNodeEqID() const;
00087 
00088     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > >::type& getWsElNodeID() const;
00089 
00091     Teuchos::ArrayRCP<double>& getCoordinates() const;
00092 
00093     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type& getCoords() const;
00094     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type& getSurfaceHeight() const;
00095     const Albany::WorksetArray<Teuchos::ArrayRCP<double> >::type& getTemperature() const;
00096     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type& getBasalFriction() const;
00097     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type& getThickness() const;
00098     const Albany::WorksetArray<Teuchos::ArrayRCP<double> >::type& getFlowFactor() const;
00099     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type& getSurfaceVelocity() const;
00100     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type& getVelocityRMS() const;
00101 
00103 
00104     void printCoords() const;
00105 
00106     Albany::StateArrays& getStateArrays() {return stateArrays;}
00107 
00109     const Albany::WorksetArray<std::string>::type&  getWsEBNames() const;
00111     const Albany::WorksetArray<int>::type&  getWsPhysIndex() const;
00112 
00113     void writeSolution(const Epetra_Vector& soln, const double time, const bool overlapped = false);
00114 
00115     Teuchos::RCP<Epetra_Vector> getSolutionField() const;
00116 
00117     int getSolutionFieldHistoryDepth() const;
00118     Teuchos::RCP<Epetra_MultiVector> getSolutionFieldHistory() const;
00119     Teuchos::RCP<Epetra_MultiVector> getSolutionFieldHistory(int maxStepCount) const;
00120     void getSolutionFieldHistory(Epetra_MultiVector &result) const;
00121 
00122     void setResidualField(const Epetra_Vector& residual);
00123 
00124     // Retrieve mesh struct
00125     Teuchos::RCP<Albany::AbstractSTKMeshStruct> getSTKMeshStruct() {return stkMeshStruct;}
00126     Teuchos::RCP<Albany::AbstractMeshStruct> getMeshStruct() const {return stkMeshStruct;}
00127 
00129     bool hasRestartSolution() const {return stkMeshStruct->hasRestartSolution();}
00130 
00132     virtual bool supportsMOR() const { return true; }
00133 
00135     double restartDataTime() const {return stkMeshStruct->restartDataTime();}
00136 
00138     void updateMesh();
00139 
00141     void transformMesh();
00142 
00144     void reNameExodusOutput(std::string& filename);
00145 
00147     int getNumDim() const { return stkMeshStruct->numDim; }
00148 
00150     int getNumEq() const { return neq; }
00151 
00153     int getOwnedDOF(const int inode, const int eq) const;
00154 
00156     int getOverlapDOF(const int inode, const int eq) const;
00157 
00159     int getGlobalDOF(const int inode, const int eq) const;
00160 
00161 
00163     // Each struct contains a latitude/longitude index and it's parametric
00164     // coordinates in an element.
00165     struct interp {
00166       std::pair<double, double> parametric_coords;
00167       std::pair<unsigned, unsigned> latitude_longitude;
00168     };
00169   private:
00170 
00172     STKDiscretization(const STKDiscretization&);
00173 
00175     STKDiscretization& operator=(const STKDiscretization&);
00176 
00177     inline int gid(const stk::mesh::Entity& node) const;
00178     inline int gid(const stk::mesh::Entity* node) const;
00179 
00180     // Copy values from STK Mesh field to given Epetra_Vector
00181     void getSolutionField(Epetra_Vector &result) const;
00182 
00183     Teuchos::RCP<Epetra_MultiVector> getSolutionFieldHistoryImpl(int stepCount) const;
00184     void getSolutionFieldHistoryImpl(Epetra_MultiVector &result) const;
00185 
00186     // Copy solution vector from Epetra_Vector into STK Mesh
00187     // Here soln is the local (non overlapped) solution
00188     void setSolutionField(const Epetra_Vector& soln);
00189 
00190     // Copy solution vector from Epetra_Vector into STK Mesh
00191     // Here soln is the local + neighbor (overlapped) solution
00192     void setOvlpSolutionField(const Epetra_Vector& soln);
00193 
00194     int nonzeroesPerRow(const int neq) const;
00195     double monotonicTimeLabel(const double time);
00196 
00198     void computeOwnedNodesAndUnknowns();
00200     void setupMLCoords();
00202     void computeOverlapNodesAndUnknowns();
00204     void computeGraphs();
00206     void computeWorksetInfo();
00208     void computeNodeSets();
00210     void computeSideSets();
00212     void setupExodusOutput();
00214     void setupNetCDFOutput();
00215     int processNetCDFOutputRequest();
00217     unsigned determine_local_side_id( const stk::mesh::Entity & elem , stk::mesh::Entity & side );
00219     Teuchos::RCP<Teuchos::FancyOStream> out;
00220 
00222     void meshToGraph();
00223 
00224     double previous_time_label;
00225 
00226   protected:
00227 
00228 
00230     stk::mesh::fem::FEMMetaData& metaData;
00231     stk::mesh::BulkData& bulkData;
00232 
00234     Teuchos::RCP<const Epetra_Comm> comm;
00235 
00237     Teuchos::RCP<Epetra_Map> node_map;
00238 
00240     Teuchos::RCP<Epetra_Map> map;
00241 
00243     Teuchos::RCP<Epetra_Map> overlap_map;
00244     Teuchos::RCP<Epetra_Map> overlap_node_map;
00245 
00247     Teuchos::RCP<Epetra_CrsGraph> graph;
00248 
00250     Teuchos::RCP<Epetra_CrsGraph> overlap_graph;
00251 
00253     unsigned int myPID;
00254 
00256     const unsigned int neq;
00257 
00259     unsigned int numMyElements;
00260 
00262     Albany::NodeSetList nodeSets;
00263     Albany::NodeSetCoordList nodeSetCoords;
00264 
00266     std::vector<Albany::SideSetList> sideSets;
00267 
00269     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type wsElNodeEqID;
00270 
00271     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > >::type wsElNodeID;
00272 
00273     mutable Teuchos::ArrayRCP<double> coordinates;
00274     Albany::WorksetArray<std::string>::type wsEBNames;
00275     Albany::WorksetArray<int>::type wsPhysIndex;
00276     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type coords;
00277     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type sHeight;
00278     Albany::WorksetArray<Teuchos::ArrayRCP<double> >::type temperature;
00279     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type basalFriction;
00280     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type thickness;
00281     Albany::WorksetArray<Teuchos::ArrayRCP<double> >::type flowFactor;
00282     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type surfaceVelocity;
00283     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type velocityRMS;
00284 
00286     WsLIDList  elemGIDws;
00287 
00288     // States: vector of length worksets of a map from field name to shards array
00289     Albany::StateArrays stateArrays;
00290 
00292     std::vector< stk::mesh::Entity * > ownednodes ;
00293     std::vector< stk::mesh::Entity * > cells ;
00294 
00296     std::vector< stk::mesh::Entity * > overlapnodes ;
00297 
00299     int numOwnedNodes;
00300     int numOverlapNodes;
00301     int numGlobalNodes;
00302 
00303     // Needed to pass coordinates to ML.
00304     Teuchos::RCP<Piro::MLRigidBodyModes> rigidBodyModes;
00305 
00306     int netCDFp;
00307     int netCDFOutputRequest;
00308     int varHeight;
00309     Albany::WorksetArray<Teuchos::ArrayRCP<std::vector<interp> > >::type interpolateData;
00310 
00311     // Storage used in periodic BCs to un-roll coordinates. Pointers saved for destructor.
00312     std::vector<double*>  toDelete;
00313 
00314     Teuchos::RCP<Albany::AbstractSTKMeshStruct> stkMeshStruct;
00315 
00316     // Used in Exodus writing capability
00317 #ifdef ALBANY_SEACAS
00318     stk::io::MeshData* mesh_data;
00319 
00320     int outputInterval;
00321 #endif
00322     bool interleavedOrdering;
00323 
00324   private:
00325 
00326     MeshGraph nodalGraph;
00327 
00328     // find the location of "value" within the first "count" locations of "vector"
00329     ssize_t in_list(const std::size_t value, std::size_t count, std::size_t *vector) {
00330 
00331       for(std::size_t i=0; i < count; i++) {
00332         if(vector[i] == value)
00333           return i;
00334       }
00335        return -1;
00336     }
00337 
00338     ssize_t in_list(const std::size_t value, std::vector<std::size_t> vector) {
00339 
00340       std::size_t count = vector.size();
00341       for(std::size_t i=0; i < count; i++) {
00342         if(vector[i] == value)
00343           return i;
00344       }
00345       return -1;
00346     }
00347 
00348     ssize_t entity_in_list(const stk::mesh::Entity *value, std::vector<stk::mesh::Entity *> vector) {
00349 
00350       std::size_t count = vector.size();
00351       for(std::size_t i=0; i < count; i++) {
00352         if(vector[i]->identifier() == value->identifier())
00353           return i;
00354       }
00355       return -1;
00356     }
00357 
00358     void printVertexConnectivity();
00359 
00360   };
00361 
00362 }
00363 
00364 #endif // ALBANY_STKDISCRETIZATION_HPP