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       20683 : CartesianMeshGenerator::validParams()
      25             : {
      26       20683 :   InputParameters params = MeshGenerator::validParams();
      27       82732 :   MooseEnum dims("1=1 2 3");
      28       82732 :   params.addRequiredParam<MooseEnum>("dim", dims, "The dimension of the mesh to be generated");
      29             : 
      30       82732 :   params.addRequiredParam<std::vector<Real>>("dx", "Intervals in the X direction");
      31       82732 :   params.addParam<std::vector<unsigned int>>(
      32             :       "ix", "Number of grids in all intervals in the X direction (default to all one)");
      33       82732 :   params.addParam<std::vector<Real>>(
      34             :       "dy", "Intervals in the Y direction (required when dim>1 otherwise ignored)");
      35       82732 :   params.addParam<std::vector<unsigned int>>(
      36             :       "iy", "Number of grids in all intervals in the Y direction (default to all one)");
      37       82732 :   params.addParam<std::vector<Real>>(
      38             :       "dz", "Intervals in the Z direction (required when dim>2 otherwise ignored)");
      39       82732 :   params.addParam<std::vector<unsigned int>>(
      40             :       "iz", "Number of grids in all intervals in the Z direction (default to all one)");
      41       82732 :   params.addParam<std::vector<unsigned int>>("subdomain_id", "Block IDs (default to all zero)");
      42       20683 :   params.addClassDescription("This CartesianMeshGenerator creates a non-uniform Cartesian mesh.");
      43       41366 :   return params;
      44       20683 : }
      45             : 
      46        3201 : CartesianMeshGenerator::CartesianMeshGenerator(const InputParameters & parameters)
      47             :   : MeshGenerator(parameters),
      48        3201 :     _dim(getParam<MooseEnum>("dim")),
      49       12804 :     _dx(getParam<std::vector<Real>>("dx"))
      50             : {
      51             :   // get all other parameters if provided and check their sizes
      52        9603 :   if (isParamValid("ix"))
      53             :   {
      54        4558 :     _ix = getParam<std::vector<unsigned int>>("ix");
      55        2279 :     if (_ix.size() != _dx.size())
      56           0 :       mooseError("ix must be in the same size of dx");
      57        7829 :     for (unsigned int i = 0; i < _ix.size(); ++i)
      58        5550 :       if (_ix[i] == 0)
      59           0 :         mooseError("ix cannot be zero");
      60             :   }
      61             :   else
      62         922 :     _ix = std::vector<unsigned int>(_dx.size(), 1);
      63             : 
      64       10133 :   for (unsigned int i = 0; i < _dx.size(); ++i)
      65        6932 :     if (_dx[i] <= 0)
      66           0 :       mooseError("dx must be greater than zero. dx: ", Moose::stringify(_dx));
      67             : 
      68        9603 :   if (isParamValid("dy"))
      69             :   {
      70        4190 :     _dy = getParam<std::vector<Real>>("dy");
      71        6480 :     for (unsigned int i = 0; i < _dy.size(); ++i)
      72        4385 :       if (_dy[i] <= 0)
      73           0 :         mooseError("dy must be greater than zero. dy: ", Moose::stringify(_dy));
      74             :   }
      75             : 
      76        9603 :   if (isParamValid("iy"))
      77             :   {
      78        3900 :     _iy = getParam<std::vector<unsigned int>>("iy");
      79        1950 :     if (_iy.size() != _dy.size())
      80           0 :       mooseError("iy must be in the same size of dy");
      81        5881 :     for (unsigned int i = 0; i < _iy.size(); ++i)
      82        3931 :       if (_iy[i] == 0)
      83           0 :         mooseError("iy cannot be zero");
      84             :   }
      85             :   else
      86        1251 :     _iy = std::vector<unsigned int>(_dy.size(), 1);
      87             : 
      88       12804 :   if (isParamValid("dz"))
      89             :   {
      90        1002 :     _dz = getParam<std::vector<Real>>("dz");
      91        2020 :     for (unsigned int i = 0; i < _dz.size(); ++i)
      92        1519 :       if (_dz[i] <= 0)
      93           0 :         mooseError("dz must be greater than zero dz: ", Moose::stringify(_dz));
      94             :   }
      95             : 
      96        9603 :   if (isParamValid("iz"))
      97             :   {
      98         904 :     _iz = getParam<std::vector<unsigned int>>("iz");
      99         452 :     if (_iz.size() != _dz.size())
     100           0 :       mooseError("iz must be in the same size of dz");
     101        1775 :     for (unsigned int i = 0; i < _iz.size(); ++i)
     102        1323 :       if (_iz[i] == 0)
     103           0 :         mooseError("iz cannot be zero");
     104             :   }
     105             :   else
     106        2749 :     _iz = std::vector<unsigned int>(_dz.size(), 1);
     107             : 
     108        9603 :   if (isParamValid("subdomain_id"))
     109             :   {
     110        3628 :     _subdomain_id = getParam<std::vector<unsigned int>>("subdomain_id");
     111        6038 :     if (isParamValid("dz") && isParamValid("dy"))
     112             :     {
     113         298 :       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        4548 :     else if (isParamValid("dy"))
     120             :     {
     121        1421 :       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          95 :       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        4161 :     if (isParamValid("dz"))
     138         203 :       _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size() * _dz.size(), 0);
     139        3552 :     else if (isParamValid("dy"))
     140         173 :       _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size(), 0);
     141             :     else
     142        1011 :       _subdomain_id = std::vector<unsigned int>(_dx.size(), 0);
     143             :   }
     144             : 
     145             :   // do dimension checks and expand block IDs for all sub-grids
     146        3201 :   switch (_dim)
     147             :   {
     148        1106 :     case 1:
     149             :     {
     150        1106 :       _nx = 0;
     151        2544 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     152        1438 :         _nx += _ix[i];
     153        1106 :       _ny = 1;
     154        1106 :       _nz = 1;
     155             : 
     156        1106 :       std::vector<unsigned int> new_id;
     157        2544 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     158        6232 :         for (unsigned int ii = 0; ii < _ix[i]; ++ii)
     159        4794 :           new_id.push_back(_subdomain_id[i]);
     160             : 
     161        1106 :       _subdomain_id = new_id;
     162             : 
     163        3318 :       if (isParamValid("dy"))
     164           0 :         mooseWarning("dy provided for 1D");
     165        3318 :       if (isParamValid("iy"))
     166           0 :         mooseWarning("iy provided for 1D");
     167        3318 :       if (isParamValid("dz"))
     168           0 :         mooseWarning("dz provided for 1D");
     169        3318 :       if (isParamValid("iz"))
     170           0 :         mooseWarning("iz provided for 1D");
     171        1106 :       break;
     172        1106 :     }
     173        1594 :     case 2:
     174             :     {
     175        1594 :       _nx = 0;
     176        5726 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     177        4132 :         _nx += _ix[i];
     178        1594 :       _ny = 0;
     179        4974 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     180        3380 :         _ny += _iy[i];
     181        1594 :       _nz = 1;
     182             : 
     183        1594 :       std::vector<unsigned int> new_id;
     184        4974 :       for (unsigned int j = 0; j < _dy.size(); ++j)
     185       11735 :         for (unsigned int jj = 0; jj < _iy[j]; ++jj)
     186       32995 :           for (unsigned int i = 0; i < _dx.size(); ++i)
     187       90413 :             for (unsigned int ii = 0; ii < _ix[i]; ++ii)
     188       65773 :               new_id.push_back(_subdomain_id[j * _dx.size() + i]);
     189             : 
     190        1594 :       _subdomain_id = new_id;
     191             : 
     192        4782 :       if (!isParamValid("dy"))
     193           0 :         mooseError("dy is not provided for 2D");
     194        4782 :       if (isParamValid("dz"))
     195           0 :         mooseWarning("dz provided for 2D");
     196        4782 :       if (isParamValid("iz"))
     197           0 :         mooseWarning("iz provided for 2D");
     198        1594 :       break;
     199        1594 :     }
     200         501 :     case 3:
     201             :     {
     202         501 :       _nx = 0;
     203        1863 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     204        1362 :         _nx += _ix[i];
     205         501 :       _ny = 0;
     206        1506 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     207        1005 :         _ny += _iy[i];
     208         501 :       _nz = 0;
     209        2020 :       for (unsigned int i = 0; i < _dz.size(); ++i)
     210        1519 :         _nz += _iz[i];
     211             : 
     212         501 :       std::vector<unsigned int> new_id;
     213        2020 :       for (unsigned int k = 0; k < _dz.size(); ++k)
     214        3618 :         for (unsigned int kk = 0; kk < _iz[k]; ++kk)
     215        6314 :           for (unsigned int j = 0; j < _dy.size(); ++j)
     216       13739 :             for (unsigned int jj = 0; jj < _iy[j]; ++jj)
     217       36680 :               for (unsigned int i = 0; i < _dx.size(); ++i)
     218       63789 :                 for (unsigned int ii = 0; ii < _ix[i]; ++ii)
     219       36633 :                   new_id.push_back(_subdomain_id[k * _dx.size() * _dy.size() + j * _dx.size() + i]);
     220             : 
     221         501 :       _subdomain_id = new_id;
     222             : 
     223        1503 :       if (!isParamValid("dy"))
     224           0 :         mooseError("dy is not provided for 3D");
     225        1503 :       if (!isParamValid("dz"))
     226           0 :         mooseError("dz is not provided for 3D");
     227         501 :       break;
     228         501 :     }
     229             :   }
     230        3201 : }
     231             : 
     232             : std::unique_ptr<MeshBase>
     233        3023 : CartesianMeshGenerator::generate()
     234             : {
     235        3023 :   auto mesh = buildMeshBaseObject();
     236             : 
     237             :   // switching on MooseEnum to generate the reference mesh
     238             :   // Note: element type is fixed
     239        3023 :   switch (_dim)
     240             :   {
     241        1077 :     case 1:
     242        1077 :       MeshTools::Generation::build_line(static_cast<UnstructuredMesh &>(*mesh), _nx, 0, _nx, EDGE2);
     243        1077 :       break;
     244        1474 :     case 2:
     245        1474 :       MeshTools::Generation::build_square(
     246        1474 :           static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, 0, _nx, 0, _ny, QUAD4);
     247        1474 :       break;
     248         472 :     case 3:
     249         472 :       MeshTools::Generation::build_cube(
     250         472 :           static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, _nz, 0, _nx, 0, _ny, 0, _nz, HEX8);
     251         472 :       break;
     252             :   }
     253             : 
     254             :   // assign block IDs
     255        6046 :   MeshBase::element_iterator el = mesh->active_elements_begin();
     256        6046 :   MeshBase::element_iterator el_end = mesh->active_elements_end();
     257      101831 :   for (; el != el_end; ++el)
     258             :   {
     259       98808 :     const Point p = (*el)->vertex_average();
     260       98808 :     unsigned int ix = std::floor(p(0));
     261       98808 :     unsigned int iy = std::floor(p(1));
     262       98808 :     unsigned int iz = std::floor(p(2));
     263       98808 :     unsigned int i = iz * _nx * _ny + iy * _nx + ix;
     264       98808 :     (*el)->subdomain_id() = _subdomain_id[i];
     265             :   }
     266             : 
     267             :   // adjust node coordinates
     268        3023 :   switch (_dim)
     269             :   {
     270        1077 :     case 1:
     271             :     {
     272             :       Real base;
     273             : 
     274        1077 :       std::vector<Real> mapx;
     275             :       // Note: the starting coordinate is zero
     276        1077 :       base = 0;
     277        1077 :       mapx.push_back(base);
     278        2476 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     279             :       {
     280        6009 :         for (unsigned int j = 1; j <= _ix[i]; ++j)
     281        4610 :           mapx.push_back(base + _dx[i] / _ix[i] * j);
     282        1399 :         base += _dx[i];
     283             :       }
     284             : 
     285        1077 :       MeshBase::node_iterator node = mesh->active_nodes_begin();
     286        1077 :       MeshBase::node_iterator node_end = mesh->active_nodes_end();
     287        6572 :       for (; node != node_end; ++node)
     288             :       {
     289        5495 :         unsigned int i = (*(*node))(0) + 0.5;
     290        5495 :         (*(*node))(0) = mapx.at(i);
     291             :       }
     292        1077 :       break;
     293        1077 :     }
     294        1474 :     case 2:
     295             :     {
     296             :       Real base;
     297             : 
     298        1474 :       std::vector<Real> mapx;
     299        1474 :       base = 0;
     300        1474 :       mapx.push_back(base);
     301        5327 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     302             :       {
     303       12069 :         for (unsigned int j = 1; j <= _ix[i]; ++j)
     304        8216 :           mapx.push_back(base + _dx[i] / _ix[i] * j);
     305        3853 :         base += _dx[i];
     306             :       }
     307             : 
     308        1474 :       std::vector<Real> mapy;
     309        1474 :       base = 0;
     310        1474 :       mapy.push_back(base);
     311        4607 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     312             :       {
     313       10873 :         for (unsigned int j = 1; j <= _iy[i]; ++j)
     314        7740 :           mapy.push_back(base + _dy[i] / _iy[i] * j);
     315        3133 :         base += _dy[i];
     316             :       }
     317             : 
     318        1474 :       MeshBase::node_iterator node = mesh->active_nodes_begin();
     319        1474 :       MeshBase::node_iterator node_end = mesh->active_nodes_end();
     320       79016 :       for (; node != node_end; ++node)
     321             :       {
     322       77542 :         unsigned int i = (*(*node))(0) + 0.5;
     323       77542 :         (*(*node))(0) = mapx.at(i);
     324       77542 :         unsigned int j = (*(*node))(1) + 0.5;
     325       77542 :         (*(*node))(1) = mapy.at(j);
     326             :       }
     327        1474 :       break;
     328        1474 :     }
     329         472 :     case 3:
     330             :     {
     331             :       Real base;
     332             : 
     333         472 :       std::vector<Real> mapx;
     334         472 :       base = 0;
     335         472 :       mapx.push_back(base);
     336        1767 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     337             :       {
     338        3069 :         for (unsigned int j = 1; j <= _ix[i]; ++j)
     339        1774 :           mapx.push_back(base + _dx[i] / _ix[i] * j);
     340        1295 :         base += _dx[i];
     341             :       }
     342             : 
     343         472 :       std::vector<Real> mapy;
     344         472 :       base = 0;
     345         472 :       mapy.push_back(base);
     346        1429 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     347             :       {
     348        3069 :         for (unsigned int j = 1; j <= _iy[i]; ++j)
     349        2112 :           mapy.push_back(base + _dy[i] / _iy[i] * j);
     350         957 :         base += _dy[i];
     351             :       }
     352             : 
     353         472 :       std::vector<Real> mapz;
     354         472 :       base = 0;
     355         472 :       mapz.push_back(base);
     356        1905 :       for (unsigned int i = 0; i < _dz.size(); ++i)
     357             :       {
     358        3415 :         for (unsigned int j = 1; j <= _iz[i]; ++j)
     359        1982 :           mapz.push_back(base + _dz[i] / _iz[i] * j);
     360        1433 :         base += _dz[i];
     361             :       }
     362             : 
     363         472 :       MeshBase::node_iterator node = mesh->active_nodes_begin();
     364         472 :       MeshBase::node_iterator node_end = mesh->active_nodes_end();
     365       65259 :       for (; node != node_end; ++node)
     366             :       {
     367       64787 :         unsigned int i = (*(*node))(0) + 0.5;
     368       64787 :         (*(*node))(0) = mapx.at(i);
     369       64787 :         unsigned int j = (*(*node))(1) + 0.5;
     370       64787 :         (*(*node))(1) = mapy.at(j);
     371       64787 :         unsigned int k = (*(*node))(2) + 0.5;
     372       64787 :         (*(*node))(2) = mapz.at(k);
     373             :       }
     374         472 :       break;
     375         472 :     }
     376             :   }
     377             : 
     378        3023 :   mesh->prepare_for_use();
     379        6046 :   return dynamic_pointer_cast<MeshBase>(mesh);
     380        3023 : }