Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "CartesianMeshGenerator.h"
      11             : #include "CastUniquePointer.h"
      12             : 
      13             : // libMesh includes
      14             : #include "libmesh/mesh_generation.h"
      15             : #include "libmesh/unstructured_mesh.h"
      16             : #include "libmesh/replicated_mesh.h"
      17             : #include "libmesh/point.h"
      18             : #include "libmesh/elem.h"
      19             : #include "libmesh/node.h"
      20             : 
      21             : registerMooseObject("MooseApp", CartesianMeshGenerator);
      22             : 
      23             : InputParameters
      24       20129 : CartesianMeshGenerator::validParams()
      25             : {
      26       20129 :   InputParameters params = MeshGenerator::validParams();
      27       20129 :   MooseEnum dims("1=1 2 3");
      28       20129 :   params.addRequiredParam<MooseEnum>("dim", dims, "The dimension of the mesh to be generated");
      29             : 
      30       20129 :   params.addRequiredParam<std::vector<Real>>("dx", "Intervals in the X direction");
      31       20129 :   params.addParam<std::vector<unsigned int>>(
      32             :       "ix", "Number of grids in all intervals in the X direction (default to all one)");
      33       20129 :   params.addParam<std::vector<Real>>(
      34             :       "dy", "Intervals in the Y direction (required when dim>1 otherwise ignored)");
      35       20129 :   params.addParam<std::vector<unsigned int>>(
      36             :       "iy", "Number of grids in all intervals in the Y direction (default to all one)");
      37       20129 :   params.addParam<std::vector<Real>>(
      38             :       "dz", "Intervals in the Z direction (required when dim>2 otherwise ignored)");
      39       20129 :   params.addParam<std::vector<unsigned int>>(
      40             :       "iz", "Number of grids in all intervals in the Z direction (default to all one)");
      41       20129 :   params.addParam<std::vector<unsigned int>>("subdomain_id", "Block IDs (default to all zero)");
      42       20129 :   params.addClassDescription("This CartesianMeshGenerator creates a non-uniform Cartesian mesh.");
      43       40258 :   return params;
      44       20129 : }
      45             : 
      46        2924 : CartesianMeshGenerator::CartesianMeshGenerator(const InputParameters & parameters)
      47             :   : MeshGenerator(parameters),
      48        2924 :     _dim(getParam<MooseEnum>("dim")),
      49        5848 :     _dx(getParam<std::vector<Real>>("dx"))
      50             : {
      51             :   // get all other parameters if provided and check their sizes
      52        2924 :   if (isParamValid("ix"))
      53             :   {
      54        2101 :     _ix = getParam<std::vector<unsigned int>>("ix");
      55        2101 :     if (_ix.size() != _dx.size())
      56           0 :       mooseError("ix must be in the same size of dx");
      57        7256 :     for (unsigned int i = 0; i < _ix.size(); ++i)
      58        5155 :       if (_ix[i] == 0)
      59           0 :         mooseError("ix cannot be zero");
      60             :   }
      61             :   else
      62         823 :     _ix = std::vector<unsigned int>(_dx.size(), 1);
      63             : 
      64        9321 :   for (unsigned int i = 0; i < _dx.size(); ++i)
      65        6397 :     if (_dx[i] <= 0)
      66           0 :       mooseError("dx must be greater than zero. dx: ", Moose::stringify(_dx));
      67             : 
      68        2924 :   if (isParamValid("dy"))
      69             :   {
      70        1932 :     _dy = getParam<std::vector<Real>>("dy");
      71        6011 :     for (unsigned int i = 0; i < _dy.size(); ++i)
      72        4079 :       if (_dy[i] <= 0)
      73           0 :         mooseError("dy must be greater than zero. dy: ", Moose::stringify(_dy));
      74             :   }
      75             : 
      76        2924 :   if (isParamValid("iy"))
      77             :   {
      78        1801 :     _iy = getParam<std::vector<unsigned int>>("iy");
      79        1801 :     if (_iy.size() != _dy.size())
      80           0 :       mooseError("iy must be in the same size of dy");
      81        5470 :     for (unsigned int i = 0; i < _iy.size(); ++i)
      82        3669 :       if (_iy[i] == 0)
      83           0 :         mooseError("iy cannot be zero");
      84             :   }
      85             :   else
      86        1123 :     _iy = std::vector<unsigned int>(_dy.size(), 1);
      87             : 
      88        2924 :   if (isParamValid("dz"))
      89             :   {
      90         448 :     _dz = getParam<std::vector<Real>>("dz");
      91        1832 :     for (unsigned int i = 0; i < _dz.size(); ++i)
      92        1384 :       if (_dz[i] <= 0)
      93           0 :         mooseError("dz must be greater than zero dz: ", Moose::stringify(_dz));
      94             :   }
      95             : 
      96        2924 :   if (isParamValid("iz"))
      97             :   {
      98         404 :     _iz = getParam<std::vector<unsigned int>>("iz");
      99         404 :     if (_iz.size() != _dz.size())
     100           0 :       mooseError("iz must be in the same size of dz");
     101        1612 :     for (unsigned int i = 0; i < _iz.size(); ++i)
     102        1208 :       if (_iz[i] == 0)
     103           0 :         mooseError("iz cannot be zero");
     104             :   }
     105             :   else
     106        2520 :     _iz = std::vector<unsigned int>(_dz.size(), 1);
     107             : 
     108        2924 :   if (isParamValid("subdomain_id"))
     109             :   {
     110        1682 :     _subdomain_id = getParam<std::vector<unsigned int>>("subdomain_id");
     111        1682 :     if (isParamValid("dz") && isParamValid("dy"))
     112             :     {
     113         272 :       if (_subdomain_id.size() != _dx.size() * _dy.size() * _dz.size())
     114           0 :         mooseError(
     115           0 :             "subdomain_id must be in the size of product of sizes of dx, dy and dz. Sizes are '" +
     116           0 :             std::to_string(_subdomain_id.size()) + "' and '" +
     117           0 :             std::to_string(_dx.size() * _dy.size() * _dz.size()) + "' respectively");
     118             :     }
     119        1410 :     else if (isParamValid("dy"))
     120             :     {
     121        1322 :       if (_subdomain_id.size() != _dx.size() * _dy.size())
     122           0 :         mooseError(
     123           0 :             "subdomain_id must be in the size of product of sizes of dx and dy. Sizes are '" +
     124           0 :             std::to_string(_subdomain_id.size()) + "' and '" +
     125           0 :             std::to_string(_dx.size() * _dy.size()) + "' respectively");
     126             :     }
     127             :     else
     128             :     {
     129          88 :       if (_subdomain_id.size() != _dx.size())
     130           0 :         mooseError("subdomain_id must be in the size of product of sizes of dx. Sizes are '" +
     131           0 :                    std::to_string(_subdomain_id.size()) + "' and '" + std::to_string(_dx.size()) +
     132             :                    "' respectively");
     133             :     }
     134             :   }
     135             :   else
     136             :   {
     137        1242 :     if (isParamValid("dz"))
     138         176 :       _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size() * _dz.size(), 0);
     139        1066 :     else if (isParamValid("dy"))
     140         162 :       _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size(), 0);
     141             :     else
     142         904 :       _subdomain_id = std::vector<unsigned int>(_dx.size(), 0);
     143             :   }
     144             : 
     145             :   // do dimension checks and expand block IDs for all sub-grids
     146        2924 :   switch (_dim)
     147             :   {
     148         992 :     case 1:
     149             :     {
     150         992 :       _nx = 0;
     151        2288 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     152        1296 :         _nx += _ix[i];
     153         992 :       _ny = 1;
     154         992 :       _nz = 1;
     155             : 
     156         992 :       std::vector<unsigned int> new_id;
     157        2288 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     158        5654 :         for (unsigned int ii = 0; ii < _ix[i]; ++ii)
     159        4358 :           new_id.push_back(_subdomain_id[i]);
     160             : 
     161         992 :       _subdomain_id = new_id;
     162             : 
     163         992 :       if (isParamValid("dy"))
     164           0 :         mooseWarning("dy provided for 1D");
     165         992 :       if (isParamValid("iy"))
     166           0 :         mooseWarning("iy provided for 1D");
     167         992 :       if (isParamValid("dz"))
     168           0 :         mooseWarning("dz provided for 1D");
     169         992 :       if (isParamValid("iz"))
     170           0 :         mooseWarning("iz provided for 1D");
     171         992 :       break;
     172         992 :     }
     173        1484 :     case 2:
     174             :     {
     175        1484 :       _nx = 0;
     176        5349 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     177        3865 :         _nx += _ix[i];
     178        1484 :       _ny = 0;
     179        4655 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     180        3171 :         _ny += _iy[i];
     181        1484 :       _nz = 1;
     182             : 
     183        1484 :       std::vector<unsigned int> new_id;
     184        4655 :       for (unsigned int j = 0; j < _dy.size(); ++j)
     185       10901 :         for (unsigned int jj = 0; jj < _iy[j]; ++jj)
     186       30830 :           for (unsigned int i = 0; i < _dx.size(); ++i)
     187       83866 :             for (unsigned int ii = 0; ii < _ix[i]; ++ii)
     188       60766 :               new_id.push_back(_subdomain_id[j * _dx.size() + i]);
     189             : 
     190        1484 :       _subdomain_id = new_id;
     191             : 
     192        1484 :       if (!isParamValid("dy"))
     193           0 :         mooseError("dy is not provided for 2D");
     194        1484 :       if (isParamValid("dz"))
     195           0 :         mooseWarning("dz provided for 2D");
     196        1484 :       if (isParamValid("iz"))
     197           0 :         mooseWarning("iz provided for 2D");
     198        1484 :       break;
     199        1484 :     }
     200         448 :     case 3:
     201             :     {
     202         448 :       _nx = 0;
     203        1684 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     204        1236 :         _nx += _ix[i];
     205         448 :       _ny = 0;
     206        1356 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     207         908 :         _ny += _iy[i];
     208         448 :       _nz = 0;
     209        1832 :       for (unsigned int i = 0; i < _dz.size(); ++i)
     210        1384 :         _nz += _iz[i];
     211             : 
     212         448 :       std::vector<unsigned int> new_id;
     213        1832 :       for (unsigned int k = 0; k < _dz.size(); ++k)
     214        3250 :         for (unsigned int kk = 0; kk < _iz[k]; ++kk)
     215        5662 :           for (unsigned int j = 0; j < _dy.size(); ++j)
     216       12474 :             for (unsigned int jj = 0; jj < _iy[j]; ++jj)
     217       33488 :               for (unsigned int i = 0; i < _dx.size(); ++i)
     218       58234 :                 for (unsigned int ii = 0; ii < _ix[i]; ++ii)
     219       33424 :                   new_id.push_back(_subdomain_id[k * _dx.size() * _dy.size() + j * _dx.size() + i]);
     220             : 
     221         448 :       _subdomain_id = new_id;
     222             : 
     223         448 :       if (!isParamValid("dy"))
     224           0 :         mooseError("dy is not provided for 3D");
     225         448 :       if (!isParamValid("dz"))
     226           0 :         mooseError("dz is not provided for 3D");
     227         448 :       break;
     228         448 :     }
     229             :   }
     230        2924 : }
     231             : 
     232             : std::unique_ptr<MeshBase>
     233        2750 : CartesianMeshGenerator::generate()
     234             : {
     235        2750 :   auto mesh = buildMeshBaseObject();
     236             : 
     237             :   // switching on MooseEnum to generate the reference mesh
     238             :   // Note: element type is fixed
     239        2750 :   switch (_dim)
     240             :   {
     241         965 :     case 1:
     242         965 :       MeshTools::Generation::build_line(static_cast<UnstructuredMesh &>(*mesh), _nx, 0, _nx, EDGE2);
     243         965 :       break;
     244        1365 :     case 2:
     245        1365 :       MeshTools::Generation::build_square(
     246        1365 :           static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, 0, _nx, 0, _ny, QUAD4);
     247        1365 :       break;
     248         420 :     case 3:
     249         420 :       MeshTools::Generation::build_cube(
     250         420 :           static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, _nz, 0, _nx, 0, _ny, 0, _nz, HEX8);
     251         420 :       break;
     252             :   }
     253             : 
     254             :   // assign block IDs
     255        5500 :   MeshBase::element_iterator el = mesh->active_elements_begin();
     256        5500 :   MeshBase::element_iterator el_end = mesh->active_elements_end();
     257       93032 :   for (; el != el_end; ++el)
     258             :   {
     259       90282 :     const Point p = (*el)->vertex_average();
     260       90282 :     unsigned int ix = std::floor(p(0));
     261       90282 :     unsigned int iy = std::floor(p(1));
     262       90282 :     unsigned int iz = std::floor(p(2));
     263       90282 :     unsigned int i = iz * _nx * _ny + iy * _nx + ix;
     264       90282 :     (*el)->subdomain_id() = _subdomain_id[i];
     265             :   }
     266             : 
     267             :   // adjust node coordinates
     268        2750 :   switch (_dim)
     269             :   {
     270         965 :     case 1:
     271             :     {
     272             :       Real base;
     273             : 
     274         965 :       std::vector<Real> mapx;
     275             :       // Note: the starting coordinate is zero
     276         965 :       base = 0;
     277         965 :       mapx.push_back(base);
     278        2224 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     279             :       {
     280        5435 :         for (unsigned int j = 1; j <= _ix[i]; ++j)
     281        4176 :           mapx.push_back(base + _dx[i] / _ix[i] * j);
     282        1259 :         base += _dx[i];
     283             :       }
     284             : 
     285         965 :       MeshBase::node_iterator node = mesh->active_nodes_begin();
     286         965 :       MeshBase::node_iterator node_end = mesh->active_nodes_end();
     287        5918 :       for (; node != node_end; ++node)
     288             :       {
     289        4953 :         unsigned int i = (*(*node))(0) + 0.5;
     290        4953 :         (*(*node))(0) = mapx.at(i);
     291             :       }
     292         965 :       break;
     293         965 :     }
     294        1365 :     case 2:
     295             :     {
     296             :       Real base;
     297             : 
     298        1365 :       std::vector<Real> mapx;
     299        1365 :       base = 0;
     300        1365 :       mapx.push_back(base);
     301        4953 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     302             :       {
     303       11178 :         for (unsigned int j = 1; j <= _ix[i]; ++j)
     304        7590 :           mapx.push_back(base + _dx[i] / _ix[i] * j);
     305        3588 :         base += _dx[i];
     306             :       }
     307             : 
     308        1365 :       std::vector<Real> mapy;
     309        1365 :       base = 0;
     310        1365 :       mapy.push_back(base);
     311        4290 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     312             :       {
     313       10050 :         for (unsigned int j = 1; j <= _iy[i]; ++j)
     314        7125 :           mapy.push_back(base + _dy[i] / _iy[i] * j);
     315        2925 :         base += _dy[i];
     316             :       }
     317             : 
     318        1365 :       MeshBase::node_iterator node = mesh->active_nodes_begin();
     319        1365 :       MeshBase::node_iterator node_end = mesh->active_nodes_end();
     320       72650 :       for (; node != node_end; ++node)
     321             :       {
     322       71285 :         unsigned int i = (*(*node))(0) + 0.5;
     323       71285 :         (*(*node))(0) = mapx.at(i);
     324       71285 :         unsigned int j = (*(*node))(1) + 0.5;
     325       71285 :         (*(*node))(1) = mapy.at(j);
     326             :       }
     327        1365 :       break;
     328        1365 :     }
     329         420 :     case 3:
     330             :     {
     331             :       Real base;
     332             : 
     333         420 :       std::vector<Real> mapx;
     334         420 :       base = 0;
     335         420 :       mapx.push_back(base);
     336        1590 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     337             :       {
     338        2757 :         for (unsigned int j = 1; j <= _ix[i]; ++j)
     339        1587 :           mapx.push_back(base + _dx[i] / _ix[i] * j);
     340        1170 :         base += _dx[i];
     341             :       }
     342             : 
     343         420 :       std::vector<Real> mapy;
     344         420 :       base = 0;
     345         420 :       mapy.push_back(base);
     346        1281 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     347             :       {
     348        2779 :         for (unsigned int j = 1; j <= _iy[i]; ++j)
     349        1918 :           mapy.push_back(base + _dy[i] / _iy[i] * j);
     350         861 :         base += _dy[i];
     351             :       }
     352             : 
     353         420 :       std::vector<Real> mapz;
     354         420 :       base = 0;
     355         420 :       mapz.push_back(base);
     356        1719 :       for (unsigned int i = 0; i < _dz.size(); ++i)
     357             :       {
     358        3053 :         for (unsigned int j = 1; j <= _iz[i]; ++j)
     359        1754 :           mapz.push_back(base + _dz[i] / _iz[i] * j);
     360        1299 :         base += _dz[i];
     361             :       }
     362             : 
     363         420 :       MeshBase::node_iterator node = mesh->active_nodes_begin();
     364         420 :       MeshBase::node_iterator node_end = mesh->active_nodes_end();
     365       58970 :       for (; node != node_end; ++node)
     366             :       {
     367       58550 :         unsigned int i = (*(*node))(0) + 0.5;
     368       58550 :         (*(*node))(0) = mapx.at(i);
     369       58550 :         unsigned int j = (*(*node))(1) + 0.5;
     370       58550 :         (*(*node))(1) = mapy.at(j);
     371       58550 :         unsigned int k = (*(*node))(2) + 0.5;
     372       58550 :         (*(*node))(2) = mapz.at(k);
     373             :       }
     374         420 :       break;
     375         420 :     }
     376             :   }
     377             : 
     378        2750 :   mesh->prepare_for_use();
     379        5500 :   return dynamic_pointer_cast<MeshBase>(mesh);
     380        2750 : }