Albany: a Trilinos-based PDE code: Topology.h 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 #if !defined(LCM_Topology_Topology_h)
00007 #define LCM_Topology_Topology_h
00008 
00009 #include <iterator>
00010 
00011 #include <stk_mesh/base/FieldData.hpp>
00012 
00013 #include "Topology_Types.h"
00014 #include "Topology_FractureCriterion.h"
00015 #include "Topology_Utils.h"
00016 
00017 namespace LCM {
00018 
00019 class Topology {
00020 public:
00029   Topology(std::string const & input_file, std::string const & output_file);
00030 
00038   Topology(RCP<Albany::AbstractDiscretization> & discretization);
00039 
00049   Topology(RCP<Albany::AbstractDiscretization> & discretization,
00050       RCP<AbstractFractureCriterion> & fracture_criterion);
00051 
00061   void
00062   setEntitiesOpen(std::map<EntityKey, bool> & open_entity_map);
00063 
00064   void
00065   setEntitiesOpen();
00066 
00076   void
00077   setEntitiesOpen(EntityVector const & fractured_faces,
00078       std::map<EntityKey, bool>& open_entity_map);
00079 
00089   enum OutputType {
00090     UNIDIRECTIONAL_UNILEVEL,
00091     UNDIRECTIONAL_MULTILEVEL,
00092     BIDIRECTIONAL_UNILEVEL,
00093     BIDIRECTIONAL_MULTILEVEL
00094   };
00095 
00096   void
00097   outputToGraphviz(
00098       std::string const & output_filename,
00099       OutputType const output_type = UNIDIRECTIONAL_UNILEVEL);
00100 
00113   void
00114   graphInitialization();
00115 
00128   void
00129   removeExtraRelations();
00130 
00139   void
00140   removeNodeRelations();
00141 
00150   void
00151   removeMultiLevelRelations();
00152 
00157   std::vector<EntityVector>
00158   getElementToNodeConnectivity();
00159 
00164   void
00165   removeElementToNodeConnectivity(std::vector<EntityVector> & v);
00166 
00175   void
00176   restoreElementToNodeConnectivity();
00177 
00181   void
00182   restoreElementToNodeConnectivity(std::vector<EntityVector> & v);
00183 
00197   EntityVector
00198   getFaceNodes(Entity * entity);
00199 
00200   EntityVector
00201   getBoundaryEntityNodes(Entity const & boundary_entity);
00202 
00206   void
00207   outputBoundary();
00208 
00212   void
00213   createBoundary();
00214 
00227   EntityVector
00228   createSurfaceElementConnectivity(Entity const & face1, Entity const & face2);
00229 
00244   void
00245   createStar(
00246       Entity & entity,
00247       std::set<EntityKey> & subgraph_entities,
00248       std::set<stkEdge, EdgeLessThan> & subgraph_edges);
00249 
00262   void
00263   splitOpenFaces();
00264 
00265   void
00266   splitOpenFaces(std::map<EntityKey, bool> & open_entity_map);
00267 
00268   void
00269   splitOpenFaces(
00270       std::map<EntityKey, bool> & open_entity_map,
00271       std::vector<EntityVector>& old_connectivity,
00272       std::vector<EntityVector>& new_connectivity);
00273 
00277   void
00278   addElement(EntityRank entity_rank);
00279 
00285   void
00286   addEntities(std::vector<size_t> & requests);
00287 
00291   void
00292   removeEntity(Entity & entity);
00293 
00297   void
00298   addRelation(Entity & source_entity, Entity & target_entity,
00299       EdgeId local_relation_id);
00300 
00304   void
00305   removeRelation(Entity & source_entity, Entity & target_entity,
00306       EdgeId local_relation_id);
00307 
00312   EntityVector
00313   getEntitiesByRank(BulkData const & mesh, EntityRank entity_rank);
00314 
00318   EntityVector::size_type
00319   getNumberEntitiesByRank(BulkData const & mesh, EntityRank entity_rank);
00320 
00324   EdgeId
00325   getLocalRelationId(Entity const & source_entity,
00326       Entity const & target_entity);
00327 
00332   int
00333   getNumberLowerRankEntities(Entity const & entity);
00334 
00340   EntityVector
00341   getDirectlyConnectedEntities(
00342       Entity const & entity,
00343       EntityRank entity_rank);
00344 
00348   bool
00349   findEntityInVector(EntityVector & entities, Entity * entity);
00350 
00361   EntityVector
00362   getBoundaryEntities(Entity const & entity, EntityRank entity_rank);
00363 
00368   bool
00369   segmentIsConnected(Entity const & segment, Entity * node);
00370 
00376   EntityVector
00377   findAdjacentSegments(Entity const & segment, Entity * node);
00378 
00383   EntityVector
00384   findCellRelations(Entity const & face);
00385 
00391   EntityVector
00392   findSegmentsFromElement(Entity const & element);
00393 
00397   bool
00398   facesShareTwoPoints(Entity const & face1, Entity const & face2);
00399 
00403   EntityVector
00404   findAdjacentSegmentsFromFace(
00405       std::vector<EntityVector> const & faces_inside_element,
00406       Entity const & face,
00407       int element_number);
00408 
00412   double*
00413   getPointerOfCoordinates(Entity * entity);
00414 
00419   EntityVector
00420   getFormerElementNodes(Entity const & element,
00421       std::vector<EntityVector> const & entities);
00422 
00428   void
00429   computeBarycentricCoordinates(
00430       EntityVector const & entities,
00431       Entity * barycenter);
00432 
00436   void
00437   barycentricSubdivision();
00438 
00442   std::vector<Entity*>
00443   getClosestNodes(std::vector<std::vector<double> > points);
00444 
00450   std::vector<Entity*>
00451   getClosestNodesOnSurface(std::vector<std::vector<double> > points);
00452 
00456   double
00457   getDistanceNodeAndPoint(Entity* node, std::vector<double> point);
00458 
00463   std::vector<std::vector<double> >
00464   getCoordinatesOfTriangle(std::vector<double> const normalToPlane);
00465 
00469   double
00470   randomNumber(double valMin, double valMax);
00471 
00475   double
00476   getDistanceBetweenNodes(Entity * node1, Entity * node2);
00477 
00482   std::vector<double>
00483   getCoordinatesOfMaxAndMin();
00484 
00489   std::vector<Entity*>
00490   MeshEdgesShortestPath();
00491 
00496   std::vector<std::vector<int> >
00497   shortestpathOnBoundaryFaces(
00498       std::vector<Entity*> const & nodes,
00499       std::vector<Entity*> const & MeshEdgesShortestPath);
00500 
00504   std::vector<std::vector<int> >
00505   shortestpath(std::vector<Entity*> const & nodes);
00506 
00510   std::vector<std::vector<int> >
00511   edgesDirections();
00512 
00516   std::vector<std::vector<int> >
00517   edgesDirectionsOuterSurface();
00518 
00519 
00523   std::vector<std::vector<int> >
00524   facesDirections();
00525 
00529   std::vector<double>
00530   facesAreas();
00531 
00536   std::vector<std::vector<int> >
00537   boundaryOperator();
00538 
00543   std::vector<std::vector<double> >
00544   outputForMpsFile();
00545 
00551   std::vector<std::vector<int> >
00552   boundaryVector(std::vector<std::vector<int> > & shortPath);
00553 
00559   std::vector<std::vector<int> >
00560   boundaryVectorOuterSurface(std::vector<std::vector<int> > & shortPath);
00561 
00567   std::vector<Entity*>
00568   MinimumSurfaceFaces(std::vector<int> VectorFromLPSolver);
00569 
00573   int
00574   NumberOfRepetitions(std::vector<Entity*> & entities, Entity * entity);
00575 
00580   std::vector<double>
00581   findCoordinates(unsigned int nodeIdentifier);
00582 
00587   void
00588   barycentricSubdivision_();
00589 
00593   void
00594   divideSegmentsHalf();
00595 
00596   void
00597   addcentroid();
00598 
00599   void
00600   connectcentroid();
00601 
00602   void
00603   addnewfaces();
00604 
00605   void
00606   connectnewfaces();
00607 
00611   void
00612   setSTKMeshStruct(RCP<Albany::AbstractSTKMeshStruct> const & sms)
00613   {stk_mesh_struct_ = sms;}
00614 
00615   RCP<Albany::AbstractSTKMeshStruct> &
00616   getSTKMeshStruct()
00617   {return stk_mesh_struct_;}
00618 
00619   void
00620   setDiscretization(RCP<Albany::AbstractDiscretization> const & d)
00621   {discretization_ = d;}
00622 
00623   RCP<Albany::AbstractDiscretization> &
00624   getDiscretization()
00625   {return discretization_;}
00626 
00627   BulkData *
00628   getBulkData()
00629   {return stk_mesh_struct_->bulkData;}
00630 
00631   stk::mesh::fem::FEMMetaData *
00632   getMetaData()
00633   {return stk_mesh_struct_->metaData;}
00634 
00635   void
00636   setCellTopology(shards::CellTopology const & ct)
00637   {cell_topology_ = ct;}
00638 
00639   shards::CellTopology &
00640   getCellTopology()
00641   {return cell_topology_;}
00642 
00643   size_t const
00644   getSpaceDimension() {return static_cast<size_t>(getSTKMeshStruct()->numDim);}
00645 
00646   EntityRank const
00647   getCellRank() {return getMetaData()->element_rank();}
00648 
00649   EntityRank const
00650   getBoundaryRank()
00651   {
00652     assert(getCellRank() > 0);
00653     return getCellRank() - 1;
00654   }
00655 
00656   IntScalarFieldType &
00657   getFractureState()
00658   {return *(stk_mesh_struct_->getFieldContainer()->getFractureState());}
00659 
00660   void
00661   setFractureCriterion(RCP<AbstractFractureCriterion> const & fc)
00662   {fracture_criterion_ = fc;}
00663 
00664   RCP<AbstractFractureCriterion> &
00665   getFractureCriterion()
00666   {return fracture_criterion_;}
00667 
00668   bool
00669   isLocalEntity(Entity const & e)
00670   {return getBulkData()->parallel_rank() == e.owner_rank();}
00671 
00672   //
00673   // Set fracture state. Do nothing for cells (elements).
00674   //
00675   void
00676   setFractureState(Entity const & e, FractureState const fs)
00677   {
00678     if (e.entity_rank() < getCellRank()) {
00679       *(stk::mesh::field_data(getFractureState(), e)) = static_cast<int>(fs);
00680     }
00681   }
00682 
00683   //
00684   // Get fracture state. Return CLOSED for cells (elements).
00685   //
00686   FractureState
00687   getFractureState(Entity const & e)
00688   {
00689     return e.entity_rank() >= getCellRank() ?
00690     CLOSED :
00691     static_cast<FractureState>(*(stk::mesh::field_data(getFractureState(), e)));
00692   }
00693 
00694   bool
00695   isInternal(Entity const & e) {
00696 
00697     assert(e.entity_rank() == getBoundaryRank());
00698 
00699     PairIterRelation
00700     relations = relations_one_up(e);
00701 
00702     size_t const
00703     number_in_edges = std::distance(relations.begin(), relations.end());
00704 
00705     assert(number_in_edges == 1 || number_in_edges == 2);
00706 
00707     return number_in_edges == 2;
00708   }
00709 
00710   bool
00711   isOpen(Entity const & e) {
00712     return getFractureState(e) == OPEN;
00713   }
00714 
00715   bool
00716   isInternalAndOpen(Entity const & e) {
00717     return isInternal(e) == true && isOpen(e) == true;
00718   }
00719 
00720   bool
00721   checkOpen(Entity const & e)
00722   {
00723     return fracture_criterion_->check(e);
00724   }
00725 
00729   void
00730   initializeFractureState();
00731 
00732 private:
00738   void
00739   createDiscretization();
00740 
00744   void
00745   setHighestIds();
00746 
00747   //
00748   //
00749   RCP<Albany::AbstractDiscretization> discretization_;
00750 
00751   RCP<Albany::AbstractSTKMeshStruct> stk_mesh_struct_;
00752 
00753   std::vector<EntityVector> connectivity_temp_;
00754 
00755   std::map<int, int> element_global_to_local_ids_;
00756 
00757   std::set<EntityPair> fractured_faces_;
00758 
00759   std::vector<int> highest_ids_;
00760 
00764   shards::CellTopology cell_topology_;
00765 
00767   RCP<AbstractFractureCriterion> fracture_criterion_;
00768 
00769 protected:
00773   Topology();
00774 };
00775 // class Topology
00776 
00777 }// namespace LCM
00778 
00779 #endif // LCM_Topology_Topology_h