CartesianMeshGenerator.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "CartesianMeshGenerator.h"
#include "CastUniquePointer.h"

// libMesh includes
#include "libmesh/mesh_generation.h"
#include "libmesh/unstructured_mesh.h"
#include "libmesh/replicated_mesh.h"
#include "libmesh/point.h"
#include "libmesh/elem.h"
#include "libmesh/node.h"

registerMooseObject("MooseApp", CartesianMeshGenerator);

InputParameters
CartesianMeshGenerator::validParams()
{
  InputParameters params = MeshGenerator::validParams();
  MooseEnum dims("1=1 2 3");
  params.addRequiredParam<MooseEnum>("dim", dims, "The dimension of the mesh to be generated");

  params.addRequiredParam<std::vector<Real>>("dx", "Intervals in the X direction");
  params.addParam<std::vector<unsigned int>>(
      "ix", "Number of grids in all intervals in the X direction (default to all one)");
  params.addParam<std::vector<Real>>(
      "max_sizes_in_x",
      "Maximum element size in the X direction. This can be specified once for all 'dx' intervals "
      "or specified for each interval. 'ix' is then determined using the ceiling of dx/size");
  params.addParam<std::vector<Real>>(
      "dy", "Intervals in the Y direction (required when dim>1 otherwise ignored)");
  params.addParam<std::vector<unsigned int>>(
      "iy", "Number of grids in all intervals in the Y direction (default to all one)");
  params.addParam<std::vector<Real>>(
      "max_sizes_in_y",
      "Maximum element size in the Y direction. This can be specified once for all 'dy' intervals "
      "or specified for each interval. 'iy' is then determined using the ceiling of dy/size");
  params.addParam<std::vector<Real>>(
      "dz", "Intervals in the Z direction (required when dim>2 otherwise ignored)");
  params.addParam<std::vector<unsigned int>>(
      "iz", "Number of grids in all intervals in the Z direction (default to all one)");
  params.addParam<std::vector<Real>>(
      "max_sizes_in_z",
      "Maximum element size in the Z direction. This can be specified once for all 'dz' intervals "
      "or specified for each interval. 'iz' is then determined using the ceiling of dz/size");
  params.addParam<std::vector<unsigned int>>("subdomain_id", "Block IDs (default to all zero)");
  params.addClassDescription("This CartesianMeshGenerator creates a non-uniform Cartesian mesh.");
  return params;
}

CartesianMeshGenerator::CartesianMeshGenerator(const InputParameters & parameters)
  : MeshGenerator(parameters),
    _dim(getParam<MooseEnum>("dim")),
    _dx(getParam<std::vector<Real>>("dx"))
{
  // get all other parameters if provided and check their sizes
  if (isParamValid("ix"))
  {
    _ix = getParam<std::vector<unsigned int>>("ix");
    if (_ix.size() != _dx.size())
      mooseError("ix must be the same size as dx");
    for (unsigned int i = 0; i < _ix.size(); ++i)
      if (_ix[i] == 0)
        mooseError("ix cannot be zero");
    if (isParamValid("max_sizes_in_x"))
      paramError("max_sizes_in_x", "Cannot be provided if 'ix' is specified");
  }
  else if (isParamValid("max_sizes_in_x"))
  {
    _ix.resize(_dx.size());
    auto max_sizes = getParam<std::vector<Real>>("max_sizes_in_x");
    if (max_sizes.size() == 1)
      max_sizes.resize(_dx.size(), max_sizes[0]);
    if (max_sizes.size() != _dx.size())
      paramError("max_sizes_in_x", "max_sizes_in_x must be size 1 or the same size as dx");
    for (const auto i : index_range(_dx))
      if (max_sizes[i] <= 0)
        paramError("max_sizes_in_x", "max_sizes_in_x must be strictly positive.");
    for (const auto i : index_range(_dx))
      _ix[i] = std::ceil(_dx[i] / max_sizes[i]);
  }
  else
    _ix = std::vector<unsigned int>(_dx.size(), 1);

  for (unsigned int i = 0; i < _dx.size(); ++i)
    if (_dx[i] <= 0)
      mooseError("dx must be greater than zero. dx: ", Moose::stringify(_dx));

  if (isParamValid("dy"))
  {
    _dy = getParam<std::vector<Real>>("dy");
    for (unsigned int i = 0; i < _dy.size(); ++i)
      if (_dy[i] <= 0)
        mooseError("dy must be greater than zero. dy: ", Moose::stringify(_dy));
  }

  if (isParamValid("iy"))
  {
    _iy = getParam<std::vector<unsigned int>>("iy");
    if (_iy.size() != _dy.size())
      mooseError("iy must be the same size as dy");
    for (unsigned int i = 0; i < _iy.size(); ++i)
      if (_iy[i] == 0)
        mooseError("iy cannot be zero");
    if (isParamValid("max_sizes_in_y"))
      paramError("max_sizes_in_y", "Cannot be provided if 'iy' is specified");
  }
  else if (isParamValid("max_sizes_in_y"))
  {
    _iy.resize(_dy.size());
    auto max_sizes = getParam<std::vector<Real>>("max_sizes_in_y");
    if (max_sizes.size() == 1)
      max_sizes.resize(_dy.size(), max_sizes[0]);
    if (max_sizes.size() != _dy.size())
      paramError("max_sizes_in_y", "max_sizes_in_y must be size 1 or the same size as dy");
    for (const auto i : index_range(_dy))
      if (max_sizes[i] <= 0)
        paramError("max_sizes_in_y", "max_sizes_in_y must be strictly positive.");
    for (const auto i : index_range(_dy))
      _iy[i] = std::ceil(_dy[i] / max_sizes[i]);
  }
  else
    _iy = std::vector<unsigned int>(_dy.size(), 1);

  if (isParamValid("dz"))
  {
    _dz = getParam<std::vector<Real>>("dz");
    for (unsigned int i = 0; i < _dz.size(); ++i)
      if (_dz[i] <= 0)
        mooseError("dz must be greater than zero dz: ", Moose::stringify(_dz));
  }

  if (isParamValid("iz"))
  {
    _iz = getParam<std::vector<unsigned int>>("iz");
    if (_iz.size() != _dz.size())
      mooseError("iz must be the same size as dz");
    for (unsigned int i = 0; i < _iz.size(); ++i)
      if (_iz[i] == 0)
        mooseError("iz cannot be zero");
    if (isParamValid("max_sizes_in_z"))
      paramError("max_sizes_in_z", "Cannot be provided if 'iz' is specified");
  }
  else if (isParamValid("max_sizes_in_z"))
  {
    _iz.resize(_dz.size());
    auto max_sizes = getParam<std::vector<Real>>("max_sizes_in_z");
    if (max_sizes.size() == 1)
      max_sizes.resize(_dz.size(), max_sizes[0]);
    if (max_sizes.size() != _dz.size())
      paramError("max_sizes_in_z", "max_sizes_in_z must be size 1 or the same size as dz");
    for (const auto i : index_range(_dz))
      if (max_sizes[i] <= 0)
        paramError("max_sizes_in_z", "max_sizes_in_z must be strictly positive.");
    for (const auto i : index_range(_dz))
      _iz[i] = std::ceil(_dz[i] / max_sizes[i]);
  }
  else
    _iz = std::vector<unsigned int>(_dz.size(), 1);

  if (isParamValid("subdomain_id"))
  {
    _subdomain_id = getParam<std::vector<unsigned int>>("subdomain_id");
    if (isParamValid("dz") && isParamValid("dy"))
    {
      if (_subdomain_id.size() != _dx.size() * _dy.size() * _dz.size())
        mooseError(
            "subdomain_id must be the size of the product of sizes of dx, dy and dz. Sizes are '" +
            std::to_string(_subdomain_id.size()) + "' and '" +
            std::to_string(_dx.size() * _dy.size() * _dz.size()) + "' respectively");
    }
    else if (isParamValid("dy"))
    {
      if (_subdomain_id.size() != _dx.size() * _dy.size())
        mooseError(
            "subdomain_id must be the size of the product of sizes of dx and dy. Sizes are '" +
            std::to_string(_subdomain_id.size()) + "' and '" +
            std::to_string(_dx.size() * _dy.size()) + "' respectively");
    }
    else
    {
      if (_subdomain_id.size() != _dx.size())
        mooseError("subdomain_id must be the size of the product of sizes of dx. Sizes are '" +
                   std::to_string(_subdomain_id.size()) + "' and '" + std::to_string(_dx.size()) +
                   "' respectively");
    }
  }
  else
  {
    if (isParamValid("dz"))
      _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size() * _dz.size(), 0);
    else if (isParamValid("dy"))
      _subdomain_id = std::vector<unsigned int>(_dx.size() * _dy.size(), 0);
    else
      _subdomain_id = std::vector<unsigned int>(_dx.size(), 0);
  }

  // do dimension checks and expand block IDs for all sub-grids
  switch (_dim)
  {
    case 1:
    {
      _nx = 0;
      for (unsigned int i = 0; i < _dx.size(); ++i)
        _nx += _ix[i];
      _ny = 1;
      _nz = 1;

      std::vector<unsigned int> new_id;
      for (unsigned int i = 0; i < _dx.size(); ++i)
        for (unsigned int ii = 0; ii < _ix[i]; ++ii)
          new_id.push_back(_subdomain_id[i]);

      _subdomain_id = new_id;

      if (isParamValid("dy"))
        mooseWarning("dy provided for 1D");
      if (isParamValid("iy"))
        mooseWarning("iy provided for 1D");
      if (isParamValid("dz"))
        mooseWarning("dz provided for 1D");
      if (isParamValid("iz"))
        mooseWarning("iz provided for 1D");
      if (isParamValid("max_sizes_in_y"))
        paramWarning("max_sizes_in_y", "Should not be provided for 1D");
      if (isParamValid("max_sizes_in_z"))
        paramWarning("max_sizes_in_z", "Should not be provided for 1D");
      break;
    }
    case 2:
    {
      _nx = 0;
      for (unsigned int i = 0; i < _dx.size(); ++i)
        _nx += _ix[i];
      _ny = 0;
      for (unsigned int i = 0; i < _dy.size(); ++i)
        _ny += _iy[i];
      _nz = 1;

      std::vector<unsigned int> new_id;
      for (unsigned int j = 0; j < _dy.size(); ++j)
        for (unsigned int jj = 0; jj < _iy[j]; ++jj)
          for (unsigned int i = 0; i < _dx.size(); ++i)
            for (unsigned int ii = 0; ii < _ix[i]; ++ii)
              new_id.push_back(_subdomain_id[j * _dx.size() + i]);

      _subdomain_id = new_id;

      if (!isParamValid("dy"))
        mooseError("dy is not provided for 2D");
      if (isParamValid("dz"))
        mooseWarning("dz provided for 2D");
      if (isParamValid("iz"))
        mooseWarning("iz provided for 2D");
      if (isParamValid("max_sizes_in_z"))
        paramWarning("max_sizes_in_z", "Should not be provided for 2D");
      break;
    }
    case 3:
    {
      _nx = 0;
      for (unsigned int i = 0; i < _dx.size(); ++i)
        _nx += _ix[i];
      _ny = 0;
      for (unsigned int i = 0; i < _dy.size(); ++i)
        _ny += _iy[i];
      _nz = 0;
      for (unsigned int i = 0; i < _dz.size(); ++i)
        _nz += _iz[i];

      std::vector<unsigned int> new_id;
      for (unsigned int k = 0; k < _dz.size(); ++k)
        for (unsigned int kk = 0; kk < _iz[k]; ++kk)
          for (unsigned int j = 0; j < _dy.size(); ++j)
            for (unsigned int jj = 0; jj < _iy[j]; ++jj)
              for (unsigned int i = 0; i < _dx.size(); ++i)
                for (unsigned int ii = 0; ii < _ix[i]; ++ii)
                  new_id.push_back(_subdomain_id[k * _dx.size() * _dy.size() + j * _dx.size() + i]);

      _subdomain_id = new_id;

      if (!isParamValid("dy"))
        mooseError("dy is not provided for 3D");
      if (!isParamValid("dz"))
        mooseError("dz is not provided for 3D");
      break;
    }
  }
}

std::unique_ptr<MeshBase>
CartesianMeshGenerator::generate()
{
  auto mesh = buildMeshBaseObject();

  // switching on MooseEnum to generate the reference mesh
  // Note: element type is fixed
  switch (_dim)
  {
    case 1:
      MeshTools::Generation::build_line(static_cast<UnstructuredMesh &>(*mesh), _nx, 0, _nx, EDGE2);
      break;
    case 2:
      MeshTools::Generation::build_square(
          static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, 0, _nx, 0, _ny, QUAD4);
      break;
    case 3:
      MeshTools::Generation::build_cube(
          static_cast<UnstructuredMesh &>(*mesh), _nx, _ny, _nz, 0, _nx, 0, _ny, 0, _nz, HEX8);
      break;
  }

  // assign block IDs
  MeshBase::element_iterator el = mesh->active_elements_begin();
  MeshBase::element_iterator el_end = mesh->active_elements_end();
  for (; el != el_end; ++el)
  {
    const Point p = (*el)->vertex_average();
    unsigned int ix = std::floor(p(0));
    unsigned int iy = std::floor(p(1));
    unsigned int iz = std::floor(p(2));
    unsigned int i = iz * _nx * _ny + iy * _nx + ix;
    (*el)->subdomain_id() = _subdomain_id[i];
  }

  // adjust node coordinates
  switch (_dim)
  {
    case 1:
    {
      Real base;

      std::vector<Real> mapx;
      // Note: the starting coordinate is zero
      base = 0;
      mapx.push_back(base);
      for (unsigned int i = 0; i < _dx.size(); ++i)
      {
        for (unsigned int j = 1; j <= _ix[i]; ++j)
          mapx.push_back(base + _dx[i] / _ix[i] * j);
        base += _dx[i];
      }

      MeshBase::node_iterator node = mesh->active_nodes_begin();
      MeshBase::node_iterator node_end = mesh->active_nodes_end();
      for (; node != node_end; ++node)
      {
        unsigned int i = (*(*node))(0) + 0.5;
        (*(*node))(0) = mapx.at(i);
      }
      break;
    }
    case 2:
    {
      Real base;

      std::vector<Real> mapx;
      base = 0;
      mapx.push_back(base);
      for (unsigned int i = 0; i < _dx.size(); ++i)
      {
        for (unsigned int j = 1; j <= _ix[i]; ++j)
          mapx.push_back(base + _dx[i] / _ix[i] * j);
        base += _dx[i];
      }

      std::vector<Real> mapy;
      base = 0;
      mapy.push_back(base);
      for (unsigned int i = 0; i < _dy.size(); ++i)
      {
        for (unsigned int j = 1; j <= _iy[i]; ++j)
          mapy.push_back(base + _dy[i] / _iy[i] * j);
        base += _dy[i];
      }

      MeshBase::node_iterator node = mesh->active_nodes_begin();
      MeshBase::node_iterator node_end = mesh->active_nodes_end();
      for (; node != node_end; ++node)
      {
        unsigned int i = (*(*node))(0) + 0.5;
        (*(*node))(0) = mapx.at(i);
        unsigned int j = (*(*node))(1) + 0.5;
        (*(*node))(1) = mapy.at(j);
      }
      break;
    }
    case 3:
    {
      Real base;

      std::vector<Real> mapx;
      base = 0;
      mapx.push_back(base);
      for (unsigned int i = 0; i < _dx.size(); ++i)
      {
        for (unsigned int j = 1; j <= _ix[i]; ++j)
          mapx.push_back(base + _dx[i] / _ix[i] * j);
        base += _dx[i];
      }

      std::vector<Real> mapy;
      base = 0;
      mapy.push_back(base);
      for (unsigned int i = 0; i < _dy.size(); ++i)
      {
        for (unsigned int j = 1; j <= _iy[i]; ++j)
          mapy.push_back(base + _dy[i] / _iy[i] * j);
        base += _dy[i];
      }

      std::vector<Real> mapz;
      base = 0;
      mapz.push_back(base);
      for (unsigned int i = 0; i < _dz.size(); ++i)
      {
        for (unsigned int j = 1; j <= _iz[i]; ++j)
          mapz.push_back(base + _dz[i] / _iz[i] * j);
        base += _dz[i];
      }

      MeshBase::node_iterator node = mesh->active_nodes_begin();
      MeshBase::node_iterator node_end = mesh->active_nodes_end();
      for (; node != node_end; ++node)
      {
        unsigned int i = (*(*node))(0) + 0.5;
        (*(*node))(0) = mapx.at(i);
        unsigned int j = (*(*node))(1) + 0.5;
        (*(*node))(1) = mapy.at(j);
        unsigned int k = (*(*node))(2) + 0.5;
        (*(*node))(2) = mapz.at(k);
      }
      break;
    }
  }

  mesh->prepare_for_use();
  return dynamic_pointer_cast<MeshBase>(mesh);
}