Cell.h

// Copyright (C) 2006-2015 Anders Logg
//
// This file is part of DOLFIN.
//
// DOLFIN is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// DOLFIN is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
// Modified by Johan Hoffman 2006.
// Modified by Andre Massing 2009.
// Modified by Garth N. Wells 2010.
// Modified by Jan Blechta 2013
// Modified by Martin Alnaes, 2015

#ifndef __CELL_H
#define __CELL_H

#include <memory>

#include "CellType.h"
#include "Mesh.h"
#include "MeshEntity.h"
#include "MeshEntityIteratorBase.h"
#include "MeshFunction.h"
#include <ufc.h>
#include <dolfin/geometry/Point.h>

namespace dolfin
{


  class Cell : public MeshEntity
  {
  public:

    Cell() : MeshEntity() {}

    Cell(const Mesh& mesh, std::size_t index)
      : MeshEntity(mesh, mesh.topology().dim(), index) {}

    ~Cell() {}

    CellType::Type type() const
    { return _mesh->type().cell_type(); }

    std::size_t num_vertices() const
    { return _mesh->type().num_vertices(); }

    std::size_t orientation() const
    { return _mesh->type().orientation(*this); }

    std::size_t orientation(const Point& up) const
    { return _mesh->type().orientation(*this, up); }

    double volume() const
    { return _mesh->type().volume(*this); }

    double h() const
    { return _mesh->type().h(*this); }

    double circumradius() const
    { return _mesh->type().circumradius(*this); }

    double inradius() const
    {
      // We would need facet areas
      _mesh->init(_mesh->type().dim() - 1);

      return _mesh->type().inradius(*this);
    }

    double radius_ratio() const
    {
      // We would need facet areas
      _mesh->init(_mesh->type().dim() - 1);

      return _mesh->type().radius_ratio(*this);
    }

    double squared_distance(const Point& point) const
    { return _mesh->type().squared_distance(*this, point); }

    double distance(const Point& point) const
    {
      return sqrt(squared_distance(point));
    }

    double normal(std::size_t facet, std::size_t i) const
    { return _mesh->type().normal(*this, facet, i); }

    Point normal(std::size_t facet) const
    { return _mesh->type().normal(*this, facet); }

    Point cell_normal() const
    { return _mesh->type().cell_normal(*this); }

    double facet_area(std::size_t facet) const
    { return _mesh->type().facet_area(*this, facet); }

    void order(const std::vector<std::int64_t>& local_to_global_vertex_indices)
    { _mesh->type().order(*this, local_to_global_vertex_indices); }

    bool ordered(const std::vector<std::int64_t>& local_to_global_vertex_indices) const
    { return _mesh->type().ordered(*this, local_to_global_vertex_indices); }

    bool contains(const Point& point) const;

    bool collides(const Point& point) const;

    bool collides(const MeshEntity& entity) const;

    std::vector<Point>
    intersection(const MeshEntity& entity) const;

    // FIXME: This function is part of a UFC transition
    void get_coordinate_dofs(std::vector<double>& coordinates) const
    {
      const MeshGeometry& geom = _mesh->geometry();
      const std::size_t gdim = geom.dim();
      const std::size_t geom_degree = geom.degree();
      const std::size_t num_vertices = this->num_vertices();
      const unsigned int* vertices = this->entities(0);

      if (geom_degree == 1)
      {
        coordinates.resize(num_vertices*gdim);
        for (std::size_t i = 0; i < num_vertices; ++i)
          for (std::size_t j = 0; j < gdim; ++j)
            coordinates[i*gdim + j] = geom.x(vertices[i])[j];
      }
      else if (geom_degree == 2)
      {
        const std::size_t tdim = _mesh->topology().dim();
        const std::size_t num_edges = this->num_entities(1);
        const unsigned int* edges = this->entities(1);

        coordinates.resize((num_vertices + num_edges)*gdim);

        for (std::size_t i = 0; i < num_vertices; ++i)
          for (std::size_t j = 0; j < gdim; j++)
            coordinates[i*gdim + j] = geom.x(vertices[i])[j];

        for (std::size_t i = 0; i < num_edges; ++i)
        {
          const std::size_t entity_index
              = (tdim == 1) ? index() : edges[i];
          const std::size_t point_index
            = geom.get_entity_index(1, 0, entity_index);
          for (std::size_t j = 0; j < gdim; ++j)
            coordinates[(i + num_vertices)*gdim + j] = geom.x(point_index)[j];
        }
      }
      else
      {
        dolfin_error("Cell.h", "get coordinate_dofs", "Unsupported mesh degree");
      }

    }

    // FIXME: This function is part of a UFC transition
    void get_vertex_coordinates(std::vector<double>& coordinates) const
    {
      const std::size_t gdim = _mesh->geometry().dim();
      const std::size_t num_vertices = this->num_vertices();
      const unsigned int* vertices = this->entities(0);
      coordinates.resize(num_vertices*gdim);
      for (std::size_t i = 0; i < num_vertices; i++)
        for (std::size_t j = 0; j < gdim; j++)
          coordinates[i*gdim + j] = _mesh->geometry().x(vertices[i])[j];
    }

    // FIXME: This function is part of a UFC transition
    void get_cell_data(ufc::cell& ufc_cell, int local_facet=-1) const
    {
      ufc_cell.geometric_dimension = _mesh->geometry().dim();
      ufc_cell.local_facet = local_facet;
      if (_mesh->cell_orientations().empty())
        ufc_cell.orientation = -1;
      else
      {
        dolfin_assert(index() < _mesh->cell_orientations().size());
        ufc_cell.orientation = _mesh->cell_orientations()[index()];
      }
      ufc_cell.mesh_identifier = mesh_id();
      ufc_cell.index = index();
    }

    // FIXME: This function is part of a UFC transition
    void get_cell_topology(ufc::cell& ufc_cell) const
    {
      const MeshTopology& topology = _mesh->topology();

      const std::size_t tdim = topology.dim();
      ufc_cell.topological_dimension = tdim;
      if (_mesh->cell_orientations().empty())
        ufc_cell.orientation = -1;
      else
      {
        dolfin_assert(index() < _mesh->cell_orientations().size());
        ufc_cell.orientation = _mesh->cell_orientations()[index()];
      }
      ufc_cell.entity_indices.resize(tdim + 1);
      for (std::size_t d = 0; d < tdim; d++)
      {
        ufc_cell.entity_indices[d].resize(num_entities(d));
        if (topology.have_global_indices(d))
        {
          const std::vector<std::int64_t>& global_indices
            = topology.global_indices(d);
          for (std::size_t i = 0; i < num_entities(d); ++i)
            ufc_cell.entity_indices[d][i] = global_indices[entities(d)[i]];
        }
        else
        {
          for (std::size_t i = 0; i < num_entities(d); ++i)
            ufc_cell.entity_indices[d][i] = entities(d)[i];
        }
      }
      ufc_cell.entity_indices[tdim].resize(1);
      if (topology.have_global_indices(tdim))
        ufc_cell.entity_indices[tdim][0] = global_index();
      else
        ufc_cell.entity_indices[tdim][0] = index();

      // FIXME: Using the local cell index is inconsistent with UFC, but
      //        necessary to make DOLFIN run
      // Local cell index
      ufc_cell.index = ufc_cell.entity_indices[tdim][0];
    }

  };

  typedef MeshEntityIteratorBase<Cell> CellIterator;

}

#endif