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       20603 : CartesianMeshGenerator::validParams()
      25             : {
      26       20603 :   InputParameters params = MeshGenerator::validParams();
      27       20603 :   MooseEnum dims("1=1 2 3");
      28       20603 :   params.addRequiredParam<MooseEnum>("dim", dims, "The dimension of the mesh to be generated");
      29             : 
      30       20603 :   params.addRequiredParam<std::vector<Real>>("dx", "Intervals in the X direction");
      31       20603 :   params.addParam<std::vector<unsigned int>>(
      32             :       "ix", "Number of grids in all intervals in the X direction (default to all one)");
      33       20603 :   params.addParam<std::vector<Real>>(
      34             :       "dy", "Intervals in the Y direction (required when dim>1 otherwise ignored)");
      35       20603 :   params.addParam<std::vector<unsigned int>>(
      36             :       "iy", "Number of grids in all intervals in the Y direction (default to all one)");
      37       20603 :   params.addParam<std::vector<Real>>(
      38             :       "dz", "Intervals in the Z direction (required when dim>2 otherwise ignored)");
      39       20603 :   params.addParam<std::vector<unsigned int>>(
      40             :       "iz", "Number of grids in all intervals in the Z direction (default to all one)");
      41       20603 :   params.addParam<std::vector<unsigned int>>("subdomain_id", "Block IDs (default to all zero)");
      42       20603 :   params.addClassDescription("This CartesianMeshGenerator creates a non-uniform Cartesian mesh.");
      43       41206 :   return params;
      44       20603 : }
      45             : 
      46        3161 : CartesianMeshGenerator::CartesianMeshGenerator(const InputParameters & parameters)
      47             :   : MeshGenerator(parameters),
      48        3161 :     _dim(getParam<MooseEnum>("dim")),
      49        6322 :     _dx(getParam<std::vector<Real>>("dx"))
      50             : {
      51             :   // get all other parameters if provided and check their sizes
      52        3161 :   if (isParamValid("ix"))
      53             :   {
      54        2265 :     _ix = getParam<std::vector<unsigned int>>("ix");
      55        2265 :     if (_ix.size() != _dx.size())
      56           0 :       mooseError("ix must be in the same size of dx");
      57        7799 :     for (unsigned int i = 0; i < _ix.size(); ++i)
      58        5534 :       if (_ix[i] == 0)
      59           0 :         mooseError("ix cannot be zero");
      60             :   }
      61             :   else
      62         896 :     _ix = std::vector<unsigned int>(_dx.size(), 1);
      63             : 
      64       10051 :   for (unsigned int i = 0; i < _dx.size(); ++i)
      65        6890 :     if (_dx[i] <= 0)
      66           0 :       mooseError("dx must be greater than zero. dx: ", Moose::stringify(_dx));
      67             : 
      68        3161 :   if (isParamValid("dy"))
      69             :   {
      70        2081 :     _dy = getParam<std::vector<Real>>("dy");
      71        6450 :     for (unsigned int i = 0; i < _dy.size(); ++i)
      72        4369 :       if (_dy[i] <= 0)
      73           0 :         mooseError("dy must be greater than zero. dy: ", Moose::stringify(_dy));
      74             :   }
      75             : 
      76        3161 :   if (isParamValid("iy"))
      77             :   {
      78        1936 :     _iy = getParam<std::vector<unsigned int>>("iy");
      79        1936 :     if (_iy.size() != _dy.size())
      80           0 :       mooseError("iy must be in the same size of dy");
      81        5851 :     for (unsigned int i = 0; i < _iy.size(); ++i)
      82        3915 :       if (_iy[i] == 0)
      83           0 :         mooseError("iy cannot be zero");
      84             :   }
      85             :   else
      86        1225 :     _iy = std::vector<unsigned int>(_dy.size(), 1);
      87             : 
      88        3161 :   if (isParamValid("dz"))
      89             :   {
      90         488 :     _dz = getParam<std::vector<Real>>("dz");
      91        1994 :     for (unsigned int i = 0; i < _dz.size(); ++i)
      92        1506 :       if (_dz[i] <= 0)
      93           0 :         mooseError("dz must be greater than zero dz: ", Moose::stringify(_dz));
      94             :   }
      95             : 
      96        3161 :   if (isParamValid("iz"))
      97             :   {
      98         439 :     _iz = getParam<std::vector<unsigned int>>("iz");
      99         439 :     if (_iz.size() != _dz.size())
     100           0 :       mooseError("iz must be in the same size of dz");
     101        1749 :     for (unsigned int i = 0; i < _iz.size(); ++i)
     102        1310 :       if (_iz[i] == 0)
     103           0 :         mooseError("iz cannot be zero");
     104             :   }
     105             :   else
     106        2722 :     _iz = std::vector<unsigned int>(_dz.size(), 1);
     107             : 
     108        3161 :   if (isParamValid("subdomain_id"))
     109             :   {
     110        1813 :     _subdomain_id = getParam<std::vector<unsigned int>>("subdomain_id");
     111        1813 :     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        1515 :     else if (isParamValid("dy"))
     120             :     {
     121        1420 :       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        1348 :     if (isParamValid("dz"))
     138         190 :       _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size() * _dz.size(), 0);
     139        1158 :     else if (isParamValid("dy"))
     140         173 :       _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size(), 0);
     141             :     else
     142         985 :       _subdomain_id = std::vector<unsigned int>(_dx.size(), 0);
     143             :   }
     144             : 
     145             :   // do dimension checks and expand block IDs for all sub-grids
     146        3161 :   switch (_dim)
     147             :   {
     148        1080 :     case 1:
     149             :     {
     150        1080 :       _nx = 0;
     151        2492 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     152        1412 :         _nx += _ix[i];
     153        1080 :       _ny = 1;
     154        1080 :       _nz = 1;
     155             : 
     156        1080 :       std::vector<unsigned int> new_id;
     157        2492 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     158        6180 :         for (unsigned int ii = 0; ii < _ix[i]; ++ii)
     159        4768 :           new_id.push_back(_subdomain_id[i]);
     160             : 
     161        1080 :       _subdomain_id = new_id;
     162             : 
     163        1080 :       if (isParamValid("dy"))
     164           0 :         mooseWarning("dy provided for 1D");
     165        1080 :       if (isParamValid("iy"))
     166           0 :         mooseWarning("iy provided for 1D");
     167        1080 :       if (isParamValid("dz"))
     168           0 :         mooseWarning("dz provided for 1D");
     169        1080 :       if (isParamValid("iz"))
     170           0 :         mooseWarning("iz provided for 1D");
     171        1080 :       break;
     172        1080 :     }
     173        1593 :     case 2:
     174             :     {
     175        1593 :       _nx = 0;
     176        5722 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     177        4129 :         _nx += _ix[i];
     178        1593 :       _ny = 0;
     179        4970 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     180        3377 :         _ny += _iy[i];
     181        1593 :       _nz = 1;
     182             : 
     183        1593 :       std::vector<unsigned int> new_id;
     184        4970 :       for (unsigned int j = 0; j < _dy.size(); ++j)
     185       11726 :         for (unsigned int jj = 0; jj < _iy[j]; ++jj)
     186       32971 :           for (unsigned int i = 0; i < _dx.size(); ++i)
     187       90359 :             for (unsigned int ii = 0; ii < _ix[i]; ++ii)
     188       65737 :               new_id.push_back(_subdomain_id[j * _dx.size() + i]);
     189             : 
     190        1593 :       _subdomain_id = new_id;
     191             : 
     192        1593 :       if (!isParamValid("dy"))
     193           0 :         mooseError("dy is not provided for 2D");
     194        1593 :       if (isParamValid("dz"))
     195           0 :         mooseWarning("dz provided for 2D");
     196        1593 :       if (isParamValid("iz"))
     197           0 :         mooseWarning("iz provided for 2D");
     198        1593 :       break;
     199        1593 :     }
     200         488 :     case 3:
     201             :     {
     202         488 :       _nx = 0;
     203        1837 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     204        1349 :         _nx += _ix[i];
     205         488 :       _ny = 0;
     206        1480 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     207         992 :         _ny += _iy[i];
     208         488 :       _nz = 0;
     209        1994 :       for (unsigned int i = 0; i < _dz.size(); ++i)
     210        1506 :         _nz += _iz[i];
     211             : 
     212         488 :       std::vector<unsigned int> new_id;
     213        1994 :       for (unsigned int k = 0; k < _dz.size(); ++k)
     214        3540 :         for (unsigned int kk = 0; kk < _iz[k]; ++kk)
     215        6184 :           for (unsigned int j = 0; j < _dy.size(); ++j)
     216       13609 :             for (unsigned int jj = 0; jj < _iy[j]; ++jj)
     217       36550 :               for (unsigned int i = 0; i < _dx.size(); ++i)
     218       63529 :                 for (unsigned int ii = 0; ii < _ix[i]; ++ii)
     219       36438 :                   new_id.push_back(_subdomain_id[k * _dx.size() * _dy.size() + j * _dx.size() + i]);
     220             : 
     221         488 :       _subdomain_id = new_id;
     222             : 
     223         488 :       if (!isParamValid("dy"))
     224           0 :         mooseError("dy is not provided for 3D");
     225         488 :       if (!isParamValid("dz"))
     226           0 :         mooseError("dz is not provided for 3D");
     227         488 :       break;
     228         488 :     }
     229             :   }
     230        3161 : }
     231             : 
     232             : std::unique_ptr<MeshBase>
     233        2986 : CartesianMeshGenerator::generate()
     234             : {
     235        2986 :   auto mesh = buildMeshBaseObject();
     236             : 
     237             :   // switching on MooseEnum to generate the reference mesh
     238             :   // Note: element type is fixed
     239        2986 :   switch (_dim)
     240             :   {
     241        1053 :     case 1:
     242        1053 :       MeshTools::Generation::build_line(static_cast<UnstructuredMesh &>(*mesh), _nx, 0, _nx, EDGE2);
     243        1053 :       break;
     244        1473 :     case 2:
     245        1473 :       MeshTools::Generation::build_square(
     246        1473 :           static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, 0, _nx, 0, _ny, QUAD4);
     247        1473 :       break;
     248         460 :     case 3:
     249         460 :       MeshTools::Generation::build_cube(
     250         460 :           static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, _nz, 0, _nx, 0, _ny, 0, _nz, HEX8);
     251         460 :       break;
     252             :   }
     253             : 
     254             :   // assign block IDs
     255        5972 :   MeshBase::element_iterator el = mesh->active_elements_begin();
     256        5972 :   MeshBase::element_iterator el_end = mesh->active_elements_end();
     257      101563 :   for (; el != el_end; ++el)
     258             :   {
     259       98577 :     const Point p = (*el)->vertex_average();
     260       98577 :     unsigned int ix = std::floor(p(0));
     261       98577 :     unsigned int iy = std::floor(p(1));
     262       98577 :     unsigned int iz = std::floor(p(2));
     263       98577 :     unsigned int i = iz * _nx * _ny + iy * _nx + ix;
     264       98577 :     (*el)->subdomain_id() = _subdomain_id[i];
     265             :   }
     266             : 
     267             :   // adjust node coordinates
     268        2986 :   switch (_dim)
     269             :   {
     270        1053 :     case 1:
     271             :     {
     272             :       Real base;
     273             : 
     274        1053 :       std::vector<Real> mapx;
     275             :       // Note: the starting coordinate is zero
     276        1053 :       base = 0;
     277        1053 :       mapx.push_back(base);
     278        2428 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     279             :       {
     280        5961 :         for (unsigned int j = 1; j <= _ix[i]; ++j)
     281        4586 :           mapx.push_back(base + _dx[i] / _ix[i] * j);
     282        1375 :         base += _dx[i];
     283             :       }
     284             : 
     285        1053 :       MeshBase::node_iterator node = mesh->active_nodes_begin();
     286        1053 :       MeshBase::node_iterator node_end = mesh->active_nodes_end();
     287        6504 :       for (; node != node_end; ++node)
     288             :       {
     289        5451 :         unsigned int i = (*(*node))(0) + 0.5;
     290        5451 :         (*(*node))(0) = mapx.at(i);
     291             :       }
     292        1053 :       break;
     293        1053 :     }
     294        1473 :     case 2:
     295             :     {
     296             :       Real base;
     297             : 
     298        1473 :       std::vector<Real> mapx;
     299        1473 :       base = 0;
     300        1473 :       mapx.push_back(base);
     301        5323 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     302             :       {
     303       12060 :         for (unsigned int j = 1; j <= _ix[i]; ++j)
     304        8210 :           mapx.push_back(base + _dx[i] / _ix[i] * j);
     305        3850 :         base += _dx[i];
     306             :       }
     307             : 
     308        1473 :       std::vector<Real> mapy;
     309        1473 :       base = 0;
     310        1473 :       mapy.push_back(base);
     311        4603 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     312             :       {
     313       10864 :         for (unsigned int j = 1; j <= _iy[i]; ++j)
     314        7734 :           mapy.push_back(base + _dy[i] / _iy[i] * j);
     315        3130 :         base += _dy[i];
     316             :       }
     317             : 
     318        1473 :       MeshBase::node_iterator node = mesh->active_nodes_begin();
     319        1473 :       MeshBase::node_iterator node_end = mesh->active_nodes_end();
     320       78966 :       for (; node != node_end; ++node)
     321             :       {
     322       77493 :         unsigned int i = (*(*node))(0) + 0.5;
     323       77493 :         (*(*node))(0) = mapx.at(i);
     324       77493 :         unsigned int j = (*(*node))(1) + 0.5;
     325       77493 :         (*(*node))(1) = mapy.at(j);
     326             :       }
     327        1473 :       break;
     328        1473 :     }
     329         460 :     case 3:
     330             :     {
     331             :       Real base;
     332             : 
     333         460 :       std::vector<Real> mapx;
     334         460 :       base = 0;
     335         460 :       mapx.push_back(base);
     336        1743 :       for (unsigned int i = 0; i < _dx.size(); ++i)
     337             :       {
     338        3021 :         for (unsigned int j = 1; j <= _ix[i]; ++j)
     339        1738 :           mapx.push_back(base + _dx[i] / _ix[i] * j);
     340        1283 :         base += _dx[i];
     341             :       }
     342             : 
     343         460 :       std::vector<Real> mapy;
     344         460 :       base = 0;
     345         460 :       mapy.push_back(base);
     346        1405 :       for (unsigned int i = 0; i < _dy.size(); ++i)
     347             :       {
     348        3045 :         for (unsigned int j = 1; j <= _iy[i]; ++j)
     349        2100 :           mapy.push_back(base + _dy[i] / _iy[i] * j);
     350         945 :         base += _dy[i];
     351             :       }
     352             : 
     353         460 :       std::vector<Real> mapz;
     354         460 :       base = 0;
     355         460 :       mapz.push_back(base);
     356        1881 :       for (unsigned int i = 0; i < _dz.size(); ++i)
     357             :       {
     358        3343 :         for (unsigned int j = 1; j <= _iz[i]; ++j)
     359        1922 :           mapz.push_back(base + _dz[i] / _iz[i] * j);
     360        1421 :         base += _dz[i];
     361             :       }
     362             : 
     363         460 :       MeshBase::node_iterator node = mesh->active_nodes_begin();
     364         460 :       MeshBase::node_iterator node_end = mesh->active_nodes_end();
     365       64689 :       for (; node != node_end; ++node)
     366             :       {
     367       64229 :         unsigned int i = (*(*node))(0) + 0.5;
     368       64229 :         (*(*node))(0) = mapx.at(i);
     369       64229 :         unsigned int j = (*(*node))(1) + 0.5;
     370       64229 :         (*(*node))(1) = mapy.at(j);
     371       64229 :         unsigned int k = (*(*node))(2) + 0.5;
     372       64229 :         (*(*node))(2) = mapz.at(k);
     373             :       }
     374         460 :       break;
     375         460 :     }
     376             :   }
     377             : 
     378        2986 :   mesh->prepare_for_use();
     379        5972 :   return dynamic_pointer_cast<MeshBase>(mesh);
     380        2986 : }