MeshGeometry.h

// Copyright (C) 2006 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 Garth N. Wells, 2008.
//
// First added:  2006-05-08
// Last changed: 2016-06-10

#ifndef __MESH_GEOMETRY_H
#define __MESH_GEOMETRY_H

#include <string>
#include <vector>
#include <dolfin/geometry/Point.h>
#include <dolfin/log/log.h>

namespace dolfin
{
  class Function;


  class MeshGeometry
  {
  public:

    MeshGeometry();

    MeshGeometry(const MeshGeometry& geometry);

    ~MeshGeometry();

    const MeshGeometry& operator= (const MeshGeometry& geometry);

    std::size_t dim() const
    { return _dim; }

    std::size_t degree() const
    { return _degree; }

    std::size_t num_vertices() const
    {
      dolfin_assert(coordinates.size() % _dim == 0);
      if (_degree > 1)
      {
        dolfin_assert(entity_offsets.size() > 1);
        dolfin_assert(entity_offsets[1].size() > 0);
        return entity_offsets[1][0];
      }
      return coordinates.size()/_dim;
    }

    std::size_t num_points() const
    {
      dolfin_assert(coordinates.size() % _dim == 0);
      return coordinates.size()/_dim;
    }

    const double* vertex_coordinates(std::size_t point_index)
    {
      dolfin_assert(point_index < num_vertices());
      return &coordinates[point_index*_dim];
    }

    const double* point_coordinates(std::size_t point_index)
    {
      dolfin_assert(point_index*_dim < coordinates.size());
      return &coordinates[point_index*_dim];
    }

    double x(std::size_t n, std::size_t i) const
    {
      dolfin_assert((n*_dim + i) < coordinates.size());
      dolfin_assert(i < _dim);
      return coordinates[n*_dim + i];
    }

    const double* x(std::size_t n) const
    {
      dolfin_assert(n*_dim < coordinates.size());
      return &coordinates[n*_dim];
    }

    std::vector<double>& x()
    { return coordinates; }

    const std::vector<double>& x() const
    { return coordinates; }

    Point point(std::size_t n) const;

    void init(std::size_t dim, std::size_t degree);

    void init_entities(const std::vector<std::size_t>& num_entities);

    std::size_t num_entity_coordinates(std::size_t entity_dim) const
    {
      // Calculate the number of points per entity for Lagrange
      // elements
      switch(entity_dim)
      {
      case 0:
        return 1;
      case 1:
        return (_degree - 1);
      case 2:
        return (_degree - 2)*(_degree - 1)/2;
      case 3:
        return (_degree - 3)*(_degree - 2)*(_degree - 1)/6;
      }
      dolfin_error("MeshGeometry.h",
                   "calculate number of points",
                   "Entity dimension out of range");
      return 0;
    }

    std::size_t get_entity_index(std::size_t entity_dim, std::size_t order,
                                 std::size_t index) const
    {
      dolfin_assert(entity_dim < entity_offsets.size());
      dolfin_assert(order < entity_offsets[entity_dim].size());
      const std::size_t idx = (entity_offsets[entity_dim][order] + index);
      dolfin_assert(idx*_dim < coordinates.size());
      return idx;
    }

    void set(std::size_t local_index, const double* x);

    std::size_t hash() const;

    std::string str(bool verbose) const;

  private:

    // Euclidean dimension
    std::size_t _dim;

    // Polynomial degree (1 = linear, 2 = quadratic etc.)
    std::size_t _degree;

    // Offsets to storage for coordinate points for each entity type
    std::vector<std::vector<std::size_t>> entity_offsets;

    // Coordinates for all points stored as a contiguous array
    std::vector<double> coordinates;

  };

}

#endif