Line data    Source code

       1             : /**********************************************************************/
       2             : /*                     DO NOT MODIFY THIS HEADER                      */
       3             : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
       4             : /*                                                                    */
       5             : /*            Copyright 2017 Battelle Energy Alliance, LLC            */
       6             : /*                        ALL RIGHTS RESERVED                         */
       7             : /**********************************************************************/
       8             : 
       9             : #include "MagpieUtils.h"
      10             : #include "MooseError.h"
      11             : 
      12             : #include "libmesh/fe_type.h"
      13             : #include "libmesh/fe_interface.h"
      14             : 
      15             : #include <sstream>
      16             : 
      17             : namespace MagpieUtils
      18             : {
      19             : 
      20             : Point
      21      718143 : randomElementPoint(const Elem & el, const Point & rnd)
      22             : {
      23      718143 :   FEType fe_type(el.default_order());
      24             : 
      25             :   // generate point on the reference element
      26      718143 :   Point ref(rnd);
      27             : 
      28             :   // we need to remap the random numbers according to the reference domain
      29             :   // as described in FEAbstract::on_reference_element
      30      718143 :   switch (el.type())
      31             :   {
      32             :     case NODEELEM:
      33             :       // those are just single points
      34           0 :       ref = 0.0;
      35           0 :       break;
      36             : 
      37             :     case EDGE2:
      38             :     case EDGE3:
      39             :     case EDGE4:
      40             :       // one dimensional elements in [-1:1]
      41           0 :       ref = Point(ref(0) * 2 - 1.0, 0.0, 0.0);
      42           0 :       break;
      43             : 
      44             :     case TRI3:
      45             :     case TRI6:
      46             :       // flip the upper triangle onto the lower triangle
      47           0 :       if (ref(0) + ref(1) > 1)
      48           0 :         ref = Point(1.0 - ref(0), 1.0 - ref(1), 0.0);
      49             :       else
      50           0 :         ref(2) = 0.0;
      51             :       break;
      52             : 
      53             :     case QUAD4:
      54             :     case QUAD8:
      55             :     case QUAD9:
      56             :       // two dimensional elements in [-1:1]^2
      57      636885 :       ref = Point(ref(0) * 2 - 1.0, ref(1) * 2 - 1.0, 0.0);
      58      636885 :       break;
      59             : 
      60             :     case TET4:
      61             :     case TET10:
      62             :       // already in the reference volume
      63           0 :       if (ref(0) + ref(1) + ref(2) < 1)
      64             :         break;
      65             : 
      66             :       // three more tets in the corners of the cube
      67           0 :       else if (ref(0) + 1 - ref(1) + 1 - ref(2) < 1)
      68           0 :         ref = Point(ref(0), 1.0 - ref(1), 1.0 - ref(2));
      69           0 :       else if (1 - ref(0) + ref(1) + 1 - ref(2) < 1)
      70           0 :         ref = Point(1.0 - ref(0), ref(1), 1.0 - ref(2));
      71           0 :       else if (1 - ref(0) + 1 - ref(1) + ref(2) < 1)
      72           0 :         ref = Point(1.0 - ref(0), 1.0 - ref(1), ref(2));
      73             : 
      74             :       // internal tet pair
      75             :       {
      76           0 :         Real s = (ref.norm_sq() - 1) / 2.0;
      77           0 :         ref -= Point(s, s, s);
      78             :       }
      79           0 :       break;
      80             : 
      81             :     case HEX8:
      82             :     case HEX20:
      83             :     case HEX27:
      84             :       // three dimensional elements in [-1:1]^3
      85       81258 :       ref = Point(ref(0) * 2 - 1.0, ref(1) * 2 - 1.0, ref(2) * 2 - 1.0);
      86       81258 :       break;
      87             : 
      88             :     case PRISM6:
      89             :     case PRISM15:
      90             :     case PRISM18:
      91             :       // flip the upper triangle onto the lower triangle
      92           0 :       if (ref(0) + ref(1) > 1)
      93           0 :         ref = Point(1.0 - ref(0), 1.0 - ref(1), ref(2) * 2 - 1.0);
      94             :       else
      95           0 :         ref(2) = ref(2) * 2 - 1.0;
      96             :       break;
      97             : 
      98           0 :     case PYRAMID5:
      99             :     case PYRAMID13:
     100             :     case PYRAMID14:
     101           0 :       mooseError("Random points on elements does not support pyramids");
     102             : 
     103           0 :     default:
     104           0 :       mooseError("Random points on elements does not support infinite elements");
     105             :   }
     106             : 
     107      718143 :   return FEInterface::map(el.dim(), fe_type, &el, ref);
     108             : }
     109             : 
     110             : /**
     111             :  * ZAID is an identifier of an isotope of the form
     112             :  * ZZAAAm, where m is the state.
     113             :  */
     114             : unsigned int
     115     1200174 : getZFromZAID(unsigned int zaid)
     116             : {
     117             :   mooseAssert(zaid <= 999999 && zaid > 9999, "ZAID " << zaid << " is invalid.");
     118     1200174 :   return zaid / 10000;
     119             : }
     120             : 
     121             : unsigned int
     122     1200174 : getAFromZAID(unsigned int zaid)
     123             : {
     124             :   mooseAssert(zaid <= 999999 && zaid > 9999, "ZAID " << zaid << " is invalid.");
     125     1200174 :   return (zaid % 10000) / 10;
     126             : }
     127             : 
     128             : const std::string &
     129           9 : neutronEnergyName(unsigned int i)
     130             : {
     131          14 :   const static std::vector<std::string> names({"Thermal", "Epithermal", "Fast", "High"});
     132           9 :   return names[i];
     133           2 : }
     134             : 
     135             : NeutronEnergyType
     136      500000 : determineNeutronType(Real energy)
     137             : {
     138      500000 :   if (energy < 0.5)
     139             :     return Thermal;
     140           0 :   else if (energy >= 0.5 && energy < 0.75e6)
     141             :     return Epithermal;
     142           0 :   else if (energy >= 0.75e6 && energy < 7.0e6)
     143           0 :     return Fast;
     144             :   return High;
     145             : }
     146             : 
     147             : } // namespace MagpieUtils