gridconstruction.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set ts=4 sw=2 et sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_UTILITIES_GRIDCONSTRUCTION_HH
#define DUNE_FUFEM_UTILITIES_GRIDCONSTRUCTION_HH
 
#include <memory>
 
#include <dune/common/fvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/grid/common/gridfactory.hh>
 
#include <dune/geometry/type.hh>
 
#include <dune/fufem/makering.hh>
#include <dune/fufem/makesphere.hh>
#include <dune/fufem/makehalfcircle.hh>
 
 
namespace Dune {
 
template <class GridType, int dim>
class GridConstruction {};
 
template <class GridType>
struct GridConstruction<GridType,1> {
 
  const static int dim=GridType::dimension;
  using ctype = typename GridType::ctype;
  using FV = FieldVector<ctype,dim>;
 
  static std::unique_ptr<GridType> createGrid(const ParameterTree& config) {
    auto type = config.get<std::string>("type");
    if (type == "interval")
      return createInterval(config);
    else
      DUNE_THROW(Exception,"There is no factory yet for dimension and type "+type);
  }
 
  static std::unique_ptr<GridType> createInterval(const ParameterTree& config) {
 
    FV lower = config.get<FV>("lowerCorner");
    FV upper = config.get<FV>("upperCorner");
 
    std::array<unsigned,dim> nElements = config.get<std::array<unsigned, dim> >("elements");
    FV size;
    for (int i=0; i<dim; i++)
      size[i] = (upper[i]-lower[i])/((ctype) nElements[i]);
 
    Dune::GridFactory<GridType> factory;
 
    // insert vertices
    for(ctype x=lower[0]; x<=upper[0]; x+=size[0])
      factory.insertVertex(FV({x}));
 
    using E = std::vector<unsigned>;
 
    for (unsigned i=0;i<nElements[0]; i++) {
      factory.insertElement(Dune::GeometryTypes::line, E({i, i+1}));
    }
 
    return std::unique_ptr<GridType>(factory.createGrid());
  }
};
 
template <class GridType>
struct GridConstruction<GridType,2> {
 
  const static int dim=GridType::dimension;
  using ctype = typename GridType::ctype;
  using FV = FieldVector<ctype,dim>;
 
public:
 
  static std::unique_ptr<GridType> createGrid(const ParameterTree& config) {
 
    auto type = config.get<std::string>("type");
    if (type == "rectangle")
      return createRectangle(config);
    else if (type == "circle")
      return createCircle(config);
    else
      DUNE_THROW(Exception,"There is no factory yet for dimension and type "+type);
  }
 
  static std::unique_ptr<GridType> createRectangle(const ParameterTree& config) {
 
    FV lower = config.get<FV>("lowerCorner");
    FV upper = config.get<FV>("upperCorner");
    std::array<unsigned,dim> nElements = config.get<std::array<unsigned, dim> >("elements");
    FV size;
    for (int i=0; i<dim; i++)
      size[i] = (upper[i]-lower[i])/((ctype) nElements[i]);
 
    Dune::GridFactory<GridType> factory;
 
    // insert vertices
    for(unsigned i=0; i<=nElements[0]; i++)
      for(unsigned j=0; j<=nElements[1]; j++)
          factory.insertVertex(FV({lower[0]+i*size[0], lower[1]+j*size[1]}));
 
    Dune::GeometryType gt = Dune::GeometryTypes::cube(dim);
    bool tetra = config.get<bool>("tetrahedral");
    if (tetra)
      gt = Dune::GeometryTypes::simplex(dim);
 
    using E = std::vector<unsigned>;
 
    // nodes in each slice
    unsigned nY = nElements[1]+1;
 
    for (unsigned i=0;i<nElements[0]; i++) {
      for (unsigned j=0;j<nElements[1]; j++) {
 
        if (tetra) {
          // add two triangles
          factory.insertElement(gt,E({i*nY+j, (i+1)*nY + j, i*nY + j+1}));
          factory.insertElement(gt,E({(i+1)*nY+j+1, i*nY + j+1, (i+1)*nY + j}));
        } else // add hexahedra
          factory.insertElement(gt,E({i*nY + j, (i+1)*nY + j,
                                      i*nY + j+1, (i+1)*nY + j+1}));
      }
    }
 
    return std::unique_ptr<GridType>(factory.createGrid());
  }
 
  static std::unique_ptr<GridType> createCircle(const ParameterTree& config) {
    return ::createCircle<GridType>(config.get<FV>("center"), config.get<ctype>("radius"));
  }
};
 
template <class GridType>
struct GridConstruction<GridType,3> {
 
  const static int dim=GridType::dimension;
  using ctype = typename GridType::ctype;
  using FV = Dune::FieldVector<ctype,dim>;
 
public:
 
  static std::unique_ptr<GridType> createGrid(const ParameterTree& config) {
 
    auto type = config.get<std::string>("type");
    if (type == "tubeSegment")
      return createTubeSegment(config);
    else if (type == "sphere")
      return createSphere(config);
    else if (type == "cuboid")
      return createCuboid(config);
    else
      DUNE_THROW(Exception,"There is no factory for dimension and type "+type);
  }
 
  static std::unique_ptr<GridType> createTubeSegment(const ParameterTree& config) {
 
    return makeTubeSegment3D<GridType>(config.get<FV>("center"), config.get<ctype>("thickness"),
                      config.get<ctype>("length"), config.get<ctype>("innerRadius"),
                      config.get<ctype>("fromAngle"), config.get<ctype>("toAngle"),
                      config.get<unsigned>("nElementRing"), config.get<unsigned>("nElementLength"),
                      config.get<bool>("closeTube"), config.get<size_t>("axis"),
                      config.get<bool>("tetrahedra"), config.get<bool>("parameterizedBoundary"));
  }
 
  static std::unique_ptr<GridType> createSphere(const ParameterTree& config) {
    return makeSphereOnOctahedron<GridType>(config.get<FV>("center"),
                                            config.get<ctype>("radius"));
  }
 
  static std::unique_ptr<GridType> createCuboid(const ParameterTree& config) {
 
    FV lower = config.get<FV>("lowerCorner");
    FV upper = config.get<FV>("upperCorner");
    std::array<unsigned,dim> nElements = config.get<std::array<unsigned, dim> >("elements");
    FV size;
    for (int i=0; i<dim; i++)
      size[i] = (upper[i]-lower[i])/((ctype) nElements[i]);
 
    Dune::GridFactory<GridType> factory;
 
    // insert vertices
    for(ctype x=lower[0]; x<=upper[0]; x+=size[0])
      for(ctype y=lower[1]; y<=upper[1]; y+=size[1])
        for(ctype z=lower[2]; z<=upper[2]; z+=size[2])
          factory.insertVertex(FV({x, y, z}));
 
    Dune::GeometryType gt;
    bool tetra = config.get<bool>("tetrahedral");
    if (tetra)
      gt = Dune::GeometryTypes::tetrahedron;
    else
      gt = Dune::GeometryTypes::cube(dim);
 
    using E = std::vector<unsigned>;
 
    // nodes in each z-y slice
    unsigned nZY = (nElements[1]+1)*(nElements[2]+1);
 
    for (unsigned i=0;i<nElements[0]; i++) {
      for (unsigned j=0;j<nElements[1]; j++) {
        for (unsigned k=0;k<nElements[2]; k++) {
 
          // collect indices of the hexahedral
          E hexa = {i*nZY + j*(nElements[2] + 1) + k, (i+1)*nZY + j*(nElements[2]+1) + k,
                    i*nZY + (j+1)*(nElements[2]+1) + k, (i+1)*nZY + (j+1)*(nElements[2]+1) + k,
                    i*nZY + j*(nElements[2]+1) + k + 1, (i+1)*nZY + j*(nElements[2]+1) + k + 1,
                    i*nZY + (j+1)*(nElements[2]+1) + k + 1, (i+1)*nZY + (j+1)*(nElements[2]+1) + k + 1};
 
          if (tetra) {
            // subdivide hexahedral into tetrahedra
            E tet(4);
 
            tet={hexa[0],hexa[1],hexa[3],hexa[7]};
            factory.insertElement(gt,tet);
            tet={hexa[0],hexa[5],hexa[1],hexa[7]};
            factory.insertElement(gt,tet);
            tet={hexa[0],hexa[4],hexa[5],hexa[7]};
            factory.insertElement(gt,tet);
            tet={hexa[0],hexa[4],hexa[6],hexa[7]};
            factory.insertElement(gt,tet);
            tet={hexa[0],hexa[6],hexa[2],hexa[7]};
            factory.insertElement(gt,tet);
            tet={hexa[0],hexa[3],hexa[2],hexa[7]};
            factory.insertElement(gt,tet);
          } else // add hexahedra
            factory.insertElement(gt,hexa);
        }
      }
    }
    return std::unique_ptr<GridType>(factory.createGrid());
  }
};
 
} /* namespace Dune */
 
#endif