MoabUserObject.h

#pragma once
 
// MOOSE includes
#include "UserObject.h"
#include "MaterialBase.h"
 
// MOAB includes
#include "moab/Core.hpp"
#include "moab/Skinner.hpp"
#include "moab/GeomTopoTool.hpp"
#include "MBTagConventions.hpp"
 
// Libmesh includes
#include <libmesh/elem.h>
#include <libmesh/enum_io_package.h>
#include <libmesh/enum_order.h>
#include <libmesh/enum_fe_family.h>
#include <libmesh/equation_systems.h>
#include <libmesh/system.h>
#include <libmesh/mesh_tools.h>
#include <libmesh/mesh_function.h>
 
struct MOABMaterialProperties{
  double rel_density;
  double temp;
};
 
 
// Forward Declarations
class MoabUserObject;
 
template <>
InputParameters validParams<MoabUserObject>();
 
class MoabUserObject : public UserObject
{
 public:
 
  MoabUserObject(const InputParameters & parameters);
 
  virtual void execute(){};
  virtual void initialize(){};
  virtual void finalize(){};
  virtual void threadJoin(const UserObject & /*uo*/){};
 
  void setProblem(FEProblemBase * problem) {
    if(problem==nullptr)
      throw std::logic_error("Problem passed is null");
    _problem_ptr = problem;
  } ;
 
  bool hasProblem(){ return !( _problem_ptr == nullptr ); };
 
  FEProblemBase& problem();
 
  void initMOAB();
 
  void initBinningData();
 
  void reset();
 
  bool update();
 
  bool setSolution(std::string var_now,std::vector< double > &results, double scaleFactor=1., bool isErr=false, bool normToVol=true);
 
  void getMaterialProperties(std::vector<std::string>& mat_names_out,
                             std::vector<double>& initial_densities,
                             std::vector<std::string>& tails,
                             std::vector<MOABMaterialProperties>& properties);
 
  std::shared_ptr<moab::Interface> moabPtr;
 
private:
 
  // Private types
 
  enum Sense { BACKWARDS=-1, FORWARDS=1};
 
  struct VolData{
    moab::EntityHandle vol;
    Sense sense;
  };
 
  // Private methods
 
  MeshBase& mesh();
 
  EquationSystems & systems();
 
  System& system(std::string var_now);
 
  // Helper methods to set MOAB database
 
  moab::ErrorCode createNodes(std::map<dof_id_type,moab::EntityHandle>& node_id_to_handle);
  void createElems(std::map<dof_id_type,moab::EntityHandle>& node_id_to_handle);
 
  moab::ErrorCode createTags();
 
  moab::ErrorCode createGroup(unsigned int id, std::string name,moab::EntityHandle& group_set);
 
  moab::ErrorCode createVol(unsigned int id,moab::EntityHandle& volume_set,moab::EntityHandle group_set);
 
  moab::ErrorCode createSurf(unsigned int id,moab::EntityHandle& surface_set, moab::Range& faces,  std::vector<VolData> & voldata);
 
  moab::ErrorCode createSurfaces(moab::Range& reversed, VolData& voldata, unsigned int& surf_id);
 
  moab::ErrorCode createSurfaceFromBox(const BoundingBox& box, VolData& voldata, unsigned int& surf_id, bool normalout, double factor=1.0);
 
  moab::ErrorCode createNodesFromBox(const BoundingBox& box,double factor,std::vector<moab::EntityHandle>& vert_handles);
 
  moab::ErrorCode createCornerTris(const std::vector<moab::EntityHandle> & verts,
                                   unsigned int corner, unsigned int v1,
                                   unsigned int v2 ,unsigned int v3,
                                   bool normalout, moab::Range &surface_tris);
 
  moab::ErrorCode createTri(const std::vector<moab::EntityHandle> & vertices,unsigned int v1, unsigned int v2 ,unsigned int v3, moab::Range &surface_tris);
 
  moab::ErrorCode updateSurfData(moab::EntityHandle surface_set,VolData data);
 
  moab::ErrorCode setTags(moab::EntityHandle ent,std::string name, std::string category, unsigned int id, int dim);
 
  moab::ErrorCode setTagData(moab::Tag tag, moab::EntityHandle ent, std::string data, unsigned int SIZE);
 
  moab::ErrorCode setTagData(moab::Tag tag, moab::EntityHandle ent, void* data);
 
  bool getTetSets(ElemType type, std::vector< std::vector<unsigned int> > &perms);
 
  moab::ErrorCode buildGraveyard(unsigned int & vol_id, unsigned int & surf_id);
 
  std::vector<Point> boxCoords(const BoundingBox& box, double factor);
 
  void findMaterials();
 
  void clearElemMaps();
 
  void addElem(dof_id_type id,moab::EntityHandle ent);
 
  void setSolution(unsigned int iSysNow, unsigned int iVarNow,std::vector< double > &results, double scaleFactor, bool isErr, bool normToVol);
 
  dof_id_type elem_to_soln_index(const Elem& elem,unsigned int iSysNow, unsigned int iVarNow);
 
  NumericVector<Number>& getSerialisedSolution(libMesh::System* sysPtr);
 
  void setMeshFunction(std::string var_name_in);
 
  double evalMeshFunction(std::shared_ptr<MeshFunction> meshFunctionPtr,
                          const Point& p);
 
  std::shared_ptr<MeshFunction> getMeshFunction(std::string var_name_in);
 
  bool sortElemsByResults();
 
  bool findSurfaces();
 
  void groupLocalElems(std::set<dof_id_type> elems, std::vector<moab::Range>& localElems);
 
  int getResultsBin(double value);
  inline int getResultsBinLin(double value);
  int getResultsBinLog(double value);
  inline int getRelDensityBin(double value);
 
  int getSortBin(int iVarBin, int iDenBin, int iMat,
                 int nVarBinsIn, int nDenBinsIn,int nMatsIn);
  int getSortBin(int iVarBin, int iDenBin, int iMat){
    return getSortBin(iVarBin,iDenBin,iMat,nVarBins,nDenBins,nMatBins);
  }
 
  int getMatBin(int iVarBin, int iDenBin, int nVarBinsIn, int nDenBinsIn);
  int getMatBin(int iVarBin, int iDenBin){
    return getMatBin(iVarBin,iDenBin,nVarBins,nDenBins);
  }
 
  void calcMidpoints();
  void calcMidpointsLin();
  void calcMidpointsLog();
  void calcDenMidpoints();
 
  void calcMidpointsLin(double var_min_in, double bin_width_in,int nbins_in,std::vector<double>& midpoints_in);
 
  Point elemCentroid(Elem& elem);
 
  void resetContainers();
 
  bool resetMOAB();
 
  bool findSurface(const moab::Range& region,moab::EntityHandle group, unsigned int & vol_id, unsigned int & surf_id,moab::EntityHandle& volume_set);
 
  bool write();
 
  void communicateDofSet(std::set<dof_id_type>& dofset);
 
  FEProblemBase * _problem_ptr;
 
  std::unique_ptr< moab::Skinner > skinner;
 
  std::unique_ptr< moab::GeomTopoTool > gtt;
 
  double lengthscale;
 
  double densityscale;
 
  std::map<dof_id_type,std::vector<moab::EntityHandle> > _id_to_elem_handles;
 
  moab::EntityHandle offset;
 
  // Data members relating to binning in temperature
 
  std::string var_name;
  bool binElems;
  bool logscale;
  double var_min;
  double var_max;
  double bin_width;
  int powMin;
  int powMax;
  unsigned int nVarBins;
  unsigned int nPow;
  unsigned int nMinor;
  std::vector<double> midpoints;
  std::vector<double> den_midpoints;
 
 
  std::string den_var_name;
  bool binByDensity;
  double rel_den_min;
  double rel_den_max;
  double rel_den_bw;
  unsigned int nDenBins;
  unsigned int nMatBins;
 
  std::vector<std::set<dof_id_type> > sortedElems;
 
  std::map<std::string, std::shared_ptr<MeshFunction> > meshFunctionPtrs;
 
  // N.B. For big problems this is going to be a memory bottleneck
  // TODO: We will need to come up with a better solution
  // Map is from systems index
  std::map<unsigned int, std::unique_ptr<NumericVector<Number> > > serial_solutions;
 
  // Materials data
 
  std::vector<std::string> mat_names;
  std::vector<std::string> openmc_mat_names;
  std::vector< std::set<SubdomainID> > mat_blocks;
  std::vector<double> initialDensities;
 
  moab::EntityHandle meshset;
 
  std::map<moab::EntityHandle, std::vector<VolData> > surfsToVols;
 
  // Some moab tags
 
  moab::Tag geometry_dimension_tag;
  moab::Tag id_tag;
  moab::Tag faceting_tol_tag;
  moab::Tag geometry_resabs_tag;
  moab::Tag category_tag;
  moab::Tag name_tag;
 
  const unsigned int nNodesPerTet = 4;
 
  // DagMC settings
  double faceting_tol;
  double geom_tol;
 
  double scalefactor_inner;
  double scalefactor_outer;
 
  // Settings to control the optional writing of surfaces to file.
 
  bool output_skins;
  bool output_full;
  std::string output_base;
  std::string output_base_full;
  unsigned int n_output;
  unsigned int n_period;
  unsigned int n_write;
  unsigned int n_its;
 
  PerfID _init_timer;
  PerfID _update_timer;
  PerfID _setsolution_timer;
 
};