Albany: a Trilinos-based PDE code: AlbPUMI_FMDBDiscretization.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 ALBPUMI_FMDBDISCRETIZATION_HPP
00008 #define ALBPUMI_FMDBDISCRETIZATION_HPP
00009 
00010 #include <vector>
00011 #include <fstream>
00012 
00013 #include "Teuchos_ParameterList.hpp"
00014 #include "Teuchos_VerboseObject.hpp"
00015 #include "Epetra_Comm.h"
00016 
00017 #include "AlbPUMI_AbstractPUMIDiscretization.hpp"
00018 #include "AlbPUMI_FMDBMeshStruct.hpp"
00019 #include "AlbPUMI_FMDBVtk.hpp"
00020 #include "AlbPUMI_FMDBExodus.hpp"
00021 
00022 #include "Piro_NullSpaceUtils.hpp" // has defn of struct that holds null space info for ML
00023 
00024 #include "Epetra_CrsMatrix.h"
00025 #include "Epetra_Vector.h"
00026 
00027 namespace AlbPUMI {
00028 
00029 template<class Output>
00030   class FMDBDiscretization : public AbstractPUMIDiscretization {
00031   public:
00032 
00034     FMDBDiscretization(
00035        Teuchos::RCP<AlbPUMI::FMDBMeshStruct> fmdbMeshStruct,
00036        const Teuchos::RCP<const Epetra_Comm>& comm,
00037        const Teuchos::RCP<Piro::MLRigidBodyModes>& rigidBodyModes = Teuchos::null);
00038 
00039 
00041     ~FMDBDiscretization();
00042 
00044     Teuchos::RCP<const Epetra_Map> getMap() const;
00045 
00047     Teuchos::RCP<const Epetra_Map> getOverlapMap() const;
00048 
00050     Teuchos::RCP<const Epetra_CrsGraph> getJacobianGraph() const;
00051 
00053     Teuchos::RCP<const Epetra_CrsGraph> getOverlapJacobianGraph() const;
00054 
00056     Teuchos::RCP<const Epetra_Map> getNodeMap() const;
00057 
00059     Teuchos::RCP<const Epetra_Map> getOverlapNodeMap() const;
00060 
00062     void setupMLCoords();
00063 
00065     const Albany::NodeSetList& getNodeSets() const { return nodeSets; };
00066     const Albany::NodeSetCoordList& getNodeSetCoords() const { return nodeSetCoords; };
00067 
00069     const Albany::SideSetList& getSideSets(const int workset) const { return sideSets[workset]; };
00070 
00072     Albany::WsLIDList& getElemGIDws() { return elemGIDws; };
00073 
00075     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type& getWsElNodeEqID() const;
00076 
00077     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > >::type& getWsElNodeID() const;
00078 
00080     Teuchos::ArrayRCP<double>& getCoordinates() const;
00081 
00082     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type& getCoords() const;
00083 
00084     // FIXME - Dummy FELIX accessor functions
00085     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type& getSurfaceHeight() const;
00086     const Albany::WorksetArray<Teuchos::ArrayRCP<double> >::type& getTemperature() const;
00087     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type& getBasalFriction() const;
00088     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type& getThickness() const;
00089     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type& getSurfaceVelocity() const;
00090     const Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type& getVelocityRMS() const;
00091     const Albany::WorksetArray<Teuchos::ArrayRCP<double> >::type& getFlowFactor() const;
00092 
00094     void printCoords() const;
00095     void debugMeshWriteNative(const Epetra_Vector& sol, const char* filename);
00096     void debugMeshWrite(const Epetra_Vector& sol, const char* filename);
00097 
00099     int getNumDim() const { return fmdbMeshStruct->numDim; }
00100 
00101     virtual Teuchos::RCP<const Epetra_Comm> getComm() const { return comm; }
00102 
00104     int getNumEq() const { return neq; }
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     void setResidualField(const Epetra_Vector& residual);
00118 
00119     // Retrieve mesh struct
00120     Teuchos::RCP<AlbPUMI::FMDBMeshStruct> getFMDBMeshStruct() {return fmdbMeshStruct;}
00121     Teuchos::RCP<Albany::AbstractMeshStruct> getMeshStruct() const {return fmdbMeshStruct;}
00122 
00124     bool hasRestartSolution() const {return fmdbMeshStruct->hasRestartSolution;}
00125 
00127     double restartDataTime() const {return fmdbMeshStruct->restartDataTime;}
00128 
00130     virtual bool supportsMOR() const { return false; }
00131 
00132     // After mesh modification, need to update the element connectivity and nodal coordinates
00133     void updateMesh();
00134 
00136     void getOwned_xyz(double **x, double **y, double **z, double **rbm,
00137                       int& nNodes, int numPDEs, int numScalar, int nullSpaceDim);
00138 
00139     // Function that transforms an FMDB mesh of a unit cube (for FELIX problems)
00140     // not supported in FMDB now
00141     void transformMesh(){}
00142 
00143     int getDOF(const int inode, const int eq) const
00144     {
00145       if (interleavedOrdering) return inode*neq + eq;
00146       else  return inode + numOwnedNodes*eq;
00147     }
00148 
00149     int getOwnedDOF(const int inode, const int eq) const
00150     {
00151       return getDOF(inode,eq);
00152     }
00153 
00154     int getOverlapDOF(const int inode, const int eq) const
00155     {
00156       return getDOF(inode,eq);
00157     }
00158 
00159     int getGlobalDOF(const int inode, const int eq) const
00160     {
00161       return getDOF(inode,eq);
00162     }
00163 
00164     // Copy field data from Epetra_Vector to APF
00165     void setField(
00166         const char* name,
00167         const Epetra_Vector& data,
00168         bool overlapped,
00169         int offset = 0);
00170     void setSplitFields(std::vector<std::string> names, std::vector<int> indices, 
00171         const Epetra_Vector& data, bool overlapped);
00172 
00173     // Copy field data from APF to Epetra_Vector
00174     void getField(
00175         const char* name,
00176         Epetra_Vector& data,
00177         bool overlapped,
00178         int offset = 0) const;
00179     void getSplitFields(std::vector<std::string> names, std::vector<int> indices,
00180         Epetra_Vector& data, bool overlapped) const;
00181 
00182     // Rename exodus output file when the problem is resized
00183     void reNameExodusOutput(const std::string& str){ meshOutput.setFileName(str);}
00184 
00185   private:
00186 
00188     FMDBDiscretization(const FMDBDiscretization&);
00189 
00191     FMDBDiscretization& operator=(const FMDBDiscretization&);
00192 
00193     // Copy solution vector from Epetra_Vector into FMDB Mesh
00194     // Here soln is the local (non overlapped) solution
00195     void setSolutionField(const Epetra_Vector& soln);
00196 
00197     int nonzeroesPerRow(const int neq) const;
00198     double monotonicTimeLabel(const double time);
00199 
00201     void computeOwnedNodesAndUnknowns();
00203     void computeOverlapNodesAndUnknowns();
00205     void computeGraphs();
00207     void computeWorksetInfo();
00209     void computeNodeSets();
00211     void computeSideSets();
00212 
00214     void copyQPScalarToAPF(unsigned nqp, QPData<double, 2>& state, apf::Field* f);
00215     void copyQPVectorToAPF(unsigned nqp, QPData<double, 3>& state, apf::Field* f);
00216     void copyQPTensorToAPF(unsigned nqp, QPData<double, 4>& state, apf::Field* f);
00217     void copyQPStatesToAPF();
00218     void removeQPStatesFromAPF();
00219 
00220     // ! Split Solution fields
00221     std::vector<std::string> solNames;
00222     std::vector<std::string> resNames;
00223     std::vector<int> solIndex;
00224 
00226     Teuchos::RCP<Teuchos::FancyOStream> out;
00227 
00228     double previous_time_label;
00229 
00230     // Transformation types for FELIX problems
00231     enum TRANSFORMTYPE {NONE, ISMIP_HOM_TEST_A};
00232     TRANSFORMTYPE transform_type;
00233 
00234   protected:
00235 
00237     Output meshOutput;
00238 
00240 
00242     Teuchos::RCP<const Epetra_Comm> comm;
00243 
00245     Teuchos::RCP<Epetra_Map> node_map;
00246 
00248     Teuchos::RCP<Epetra_Map> map;
00249 
00251     Teuchos::RCP<Epetra_Map> overlap_map;
00252     Teuchos::RCP<Epetra_Map> overlap_node_map;
00253 
00255     Teuchos::RCP<Epetra_CrsGraph> graph;
00256 
00258     Teuchos::RCP<Epetra_CrsGraph> overlap_graph;
00259 
00261     const unsigned int neq;
00262 
00264     Albany::NodeSetList nodeSets;
00265     Albany::NodeSetCoordList nodeSetCoords;
00266 
00268     std::vector<Albany::SideSetList> sideSets;
00269 
00271     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > > >::type wsElNodeEqID;
00272 
00273     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> > >::type wsElNodeID;
00274 
00275     mutable Teuchos::ArrayRCP<double> coordinates;
00276     Albany::WorksetArray<std::string>::type wsEBNames;
00277     Albany::WorksetArray<int>::type wsPhysIndex;
00278     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type coords;
00279 
00280     // FELIX unused variables (FIXME)
00281     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type sHeight;
00282     Albany::WorksetArray<Teuchos::ArrayRCP<double> >::type temperature;
00283     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type basalFriction;
00284     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double> > >::type thickness;
00285     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type surfaceVelocity;
00286     Albany::WorksetArray<Teuchos::ArrayRCP<Teuchos::ArrayRCP<double*> > >::type velocityRMS;
00287     Albany::WorksetArray<Teuchos::ArrayRCP<double> >::type flowFactor;
00288 
00290     Albany::WsLIDList  elemGIDws;
00291 
00292     // States: vector of length num worksets of a map from field name to shards array
00293     Albany::StateArrays stateArrays;
00294 
00296     apf::DynamicArray<apf::Node> nodes;
00297 
00299     int numOwnedNodes;
00300     int numOverlapNodes;
00301     int numGlobalNodes;
00302 
00303     // Coordinate vector in format needed by ML. Need to own memory here.
00304     double *xx, *yy, *zz, *rr;
00305     bool allocated_xyz;
00306 
00307     Teuchos::RCP<AlbPUMI::FMDBMeshStruct> fmdbMeshStruct;
00308 
00309     bool interleavedOrdering;
00310 
00311     std::vector< std::vector<apf::MeshEntity*> > buckets; // bucket of elements
00312 
00313     // storage to save the node coordinates of the nodesets visible to this PE
00314     std::map<std::string, std::vector<double> > nodeset_node_coords;
00315 
00316     // Needed to pass coordinates to ML. 
00317     Teuchos::RCP<Piro::MLRigidBodyModes> rigidBodyModes;
00318 
00319     // counter for limiting data writes to output file
00320     int outputInterval;
00321 
00322   };
00323 
00324 }
00325 
00326 // Define macro for explicit template instantiation
00327 #define FMDB_INSTANTIATE_TEMPLATE_CLASS_VTK(name) \
00328   template class name<AlbPUMI::FMDBVtk>;
00329 #define FMDB_INSTANTIATE_TEMPLATE_CLASS_EXODUS(name) \
00330   template class name<AlbPUMI::FMDBExodus>;
00331 
00332 #define FMDB_INSTANTIATE_TEMPLATE_CLASS(name) \
00333   FMDB_INSTANTIATE_TEMPLATE_CLASS_VTK(name) \
00334   FMDB_INSTANTIATE_TEMPLATE_CLASS_EXODUS(name)
00335 
00336 #endif // ALBANY_FMDBDISCRETIZATION_HPP