mcell/libmcell/api/geometry_utils.cpp

/******************************************************************************
 *
 * Copyright (C) 2020 by
 * The Salk Institute for Biological Studies
 *
 * Use of this source code is governed by an MIT-style
 * license that can be found in the LICENSE file or at
 * https://opensource.org/licenses/MIT.
 *
******************************************************************************/

#include "generated/gen_geometry_utils.h"

#include "api/geometry_object.h"
#include "api/model.h"

#include "world.h"
#include "geometry.h"
#include "partition.h"

#include <vtkPlatonicSolidSource.h>
#include <vtkPolyData.h>
#include <vtkSmartPointer.h>

using namespace std;

namespace MCell {
namespace API {

namespace geometry_utils {

std::shared_ptr<GeometryObject> create_box(
    const std::string& name,
    const double edge_dimension,
    const std::vector<double> xyz_dimensions) {

  string err_msg = S("Exactly one ") + NAME_EDGE_DIMENSION + " or individual x, y, z dimensions can be set.";
  double x, y, z;

  if (is_set(edge_dimension)) {
    if (is_set(xyz_dimensions)) {
      throw ValueError(err_msg);
    }
    x = edge_dimension / 2;
    y = edge_dimension / 2;
    z = edge_dimension / 2;
  }
  else {
    if (xyz_dimensions.size() != 3) {
      throw ValueError(err_msg);
    }

    x = xyz_dimensions[0] / 2;
    y = xyz_dimensions[1] / 2;
    z = xyz_dimensions[2] / 2;
  }

  // using the same ordering of vertices and faces as what is generated by
  // CellBlender, changing this impacts surface molecule releases

  vector<vector<double>> vertex_list {
      {-x, -y, -z},
      {-x, -y,  z},
      {-x,  y, -z},
      {-x,  y,  z},
      { x, -y, -z},
      { x, -y,  z},
      { x,  y, -z},
      { x,  y,  z}
  };

  vector<vector<int>> wall_list {
      {1, 2, 0},
      {3, 6, 2},
      {7, 4, 6},
      {5, 0, 4},
      {6, 0, 2},
      {3, 5, 7},
      {1, 3, 2},
      {3, 7, 6},
      {7, 5, 4},
      {5, 1, 0},
      {6, 4, 0},
      {3, 1, 5}
  };

  auto res = make_shared<GeometryObject>(
      name,
      vertex_list,
      wall_list
  );

  return res;
}



// used code from the link below due to absence of any reasonable library
// https://github.com/caosdoar/spheres/blob/master/src/spheres.cpp
// published under MIT license

struct Edge
{
  uint32_t v0;
  uint32_t v1;

  Edge(uint32_t v0, uint32_t v1)
    : v0(v0 < v1 ? v0 : v1)
    , v1(v0 < v1 ? v1 : v0)
  {
  }

  bool operator <(const Edge &rhs) const
  {
    return v0 < rhs.v0 || (v0 == rhs.v0 && v1 < rhs.v1);
  }
};


static Vec3 normalize(const Vec3 &a)
{
  const double lrcp = 1.0 / std::sqrt(dot(a, a));
  return Vec3(a.x * lrcp, a.y * lrcp, a.z * lrcp);
}


static void init_icosphere(
    vector<Vec3>& vertex_list,
    vector<vector<int>>& wall_list) {

  vertex_list.clear();
  wall_list.clear();

  const double t = (1.0 + std::sqrt(5.0)) / 2.0;

  // Vertices
  vertex_list.push_back(normalize(Vec3(-1.0,  t, 0.0)));
  vertex_list.push_back(normalize(Vec3( 1.0,  t, 0.0)));
  vertex_list.push_back(normalize(Vec3(-1.0, -t, 0.0)));
  vertex_list.push_back(normalize(Vec3( 1.0, -t, 0.0)));
  vertex_list.push_back(normalize(Vec3(0.0, -1.0,  t)));
  vertex_list.push_back(normalize(Vec3(0.0,  1.0,  t)));
  vertex_list.push_back(normalize(Vec3(0.0, -1.0, -t)));
  vertex_list.push_back(normalize(Vec3(0.0,  1.0, -t)));
  vertex_list.push_back(normalize(Vec3( t, 0.0, -1.0)));
  vertex_list.push_back(normalize(Vec3( t, 0.0,  1.0)));
  vertex_list.push_back(normalize(Vec3(-t, 0.0, -1.0)));
  vertex_list.push_back(normalize(Vec3(-t, 0.0,  1.0)));

  // Faces
  wall_list.push_back(vector<int>{0, 11, 5});
  wall_list.push_back(vector<int>{0, 5, 1});
  wall_list.push_back(vector<int>{0, 1, 7});
  wall_list.push_back(vector<int>{0, 7, 10});
  wall_list.push_back(vector<int>{0, 10, 11});
  wall_list.push_back(vector<int>{1, 5, 9});
  wall_list.push_back(vector<int>{5, 11, 4});
  wall_list.push_back(vector<int>{11, 10, 2});
  wall_list.push_back(vector<int>{10, 7, 6});
  wall_list.push_back(vector<int>{7, 1, 8});
  wall_list.push_back(vector<int>{3, 9, 4});
  wall_list.push_back(vector<int>{3, 4, 2});
  wall_list.push_back(vector<int>{3, 2, 6});
  wall_list.push_back(vector<int>{3, 6, 8});
  wall_list.push_back(vector<int>{3, 8, 9});
  wall_list.push_back(vector<int>{4, 9, 5});
  wall_list.push_back(vector<int>{2, 4, 11});
  wall_list.push_back(vector<int>{6, 2, 10});
  wall_list.push_back(vector<int>{8, 6, 7});
  wall_list.push_back(vector<int>{9, 8, 1});
}


int subdivide_icosphere_edge(
    int f0, int f1,
    const Vec3& v0, const Vec3& v1,
    vector<Vec3>& vertex_list_out,
    std::map<Edge, uint32_t> &io_divisions
    )
{
  // required to make the mesh watertight
  const Edge edge(f0, f1);
  auto it = io_divisions.find(edge);
  if (it != io_divisions.end()) {
    return it->second;
  }

  const Vec3 v = normalize(Vec3(0.5) * (v0 + v1));
  const uint32_t f = vertex_list_out.size();
  vertex_list_out.push_back(v);
  io_divisions[edge] = f;
  return f;
}


static void subdivide_icosphere_mesh(
    const vector<Vec3>& vertex_list_in,
    const vector<vector<int>>& wall_list_in,
    vector<Vec3>& vertex_list_out,
    vector<vector<int>>& wall_list_out) {

  vertex_list_out = vertex_list_in;

  std::map<Edge, uint32_t> divisions; // Edge -> new vertex

  for (size_t i = 0; i < wall_list_in.size(); ++i)
  {
    const int f0 = wall_list_in[i][0];
    const int f1 = wall_list_in[i][1];
    const int f2 = wall_list_in[i][2];

    const Vec3& v0 = vertex_list_in[f0];
    const Vec3& v1 = vertex_list_in[f1];
    const Vec3& v2 = vertex_list_in[f2];

    const int f3 = subdivide_icosphere_edge(f0, f1, v0, v1, vertex_list_out, divisions);
    const int f4 = subdivide_icosphere_edge(f1, f2, v1, v2, vertex_list_out, divisions);
    const int f5 = subdivide_icosphere_edge(f2, f0, v2, v0, vertex_list_out, divisions);

    wall_list_out.push_back(vector<int>{f0, f3, f5});
    wall_list_out.push_back(vector<int>{f3, f1, f4});
    wall_list_out.push_back(vector<int>{f4, f2, f5});
    wall_list_out.push_back(vector<int>{f3, f4, f5});
  }
}


std::shared_ptr<GeometryObject> create_icosphere(
    const std::string& name, const double radius, const int subdivisions) {

  if (subdivisions < 1) {
    throw ValueError(S("Value of parameter ") + NAME_SUBDIVISIONS + " must be greater or equal to 1.");
  }

  if (subdivisions > 8) {
    throw ValueError(S("Value of parameter ") + NAME_SUBDIVISIONS + " must be less or equal to 8.");
  }

  vector<Vec3> vertex_list;
  vector<vector<int>> wall_list;

  init_icosphere(vertex_list, wall_list);

  for (int i = 0; i < subdivisions - 1; i++) {
    vector<Vec3> vertex_list_new;
    vector<vector<int>> wall_list_new;
    subdivide_icosphere_mesh(vertex_list, wall_list, vertex_list_new, wall_list_new);

    vertex_list.swap(vertex_list_new);
    wall_list.swap(wall_list_new);
  }

  // convert vertices to the required format and scale
  vector<vector<double>> double_vertices;
  for (const Vec3& v: vertex_list) {
    double_vertices.push_back(vector<double>{v.x * radius, v.y * radius, v.z * radius});
  }

  auto res = make_shared<GeometryObject>(
      name,
      double_vertices,
      wall_list
  );

  return res;
}

#if 0
// keeping this code as an example on how to transform vtkPolyData
std::shared_ptr<GeometryObject> create_sphere(
    const std::string& name, const double radius, const int resolution) {

  vtkNew<vtkSphereSource> sphere;
  sphere->SetCenter(0.0, 0.0, 0.0);
  sphere->SetRadius(radius);
  sphere->SetPhiResolution(resolution);
  sphere->SetThetaResolution(resolution);
  sphere->Update();

  vtkSmartPointer<vtkPolyData> polydata = sphere->GetOutput();

  vtkPoints* points = polydata->GetPoints();

  vector<vector<double>> vertex_list;
  int num_vertices = points->GetNumberOfPoints();
  for (int i = 0; i < num_vertices; i++) {
    double pt[3];
    points->GetPoint(i, pt);
    vertex_list.push_back(vector<double>{pt[0], pt[1], pt[2]});
  }

  vector<vector<int>> wall_list;
  int num_walls = polydata->GetNumberOfCells();
  for (int i = 0; i < num_walls; i++) {
    vtkCell* cell = polydata->GetCell(i);

    vector<int> wall;
    for (int k = 0; k < 3; k++) {
      wall.push_back(cell->GetPointId(k));
    }
    wall_list.push_back(wall);
  }

  auto res = make_shared<GeometryObject>(
      name,
      vertex_list,
      wall_list
  );

  return res;
}
#endif


void validate_volumetric_mesh(
    std::shared_ptr<Model> model,
    std::shared_ptr<GeometryObject> geometry_object) {

  if (!model->is_initialized()) {
    throw RuntimeError(S("Function ") + NAME_VALIDATE_VOLUMETRIC_MESH + " may be called only after model initialization.");
  }

  if (geometry_object->geometry_object_id == GEOMETRY_OBJECT_ID_INVALID) {
    throw RuntimeError(S("Function ") + NAME_VALIDATE_VOLUMETRIC_MESH + " may be called on " + NAME_CLASS_GEOMETRY_OBJECT + " used during model initialization.");
  }

  const Partition& p = model->get_world()->get_partition(PARTITION_ID_INITIAL);

  const MCell::GeometryObject& obj = p.get_geometry_object(geometry_object->geometry_object_id);
  string msg = obj.validate_volumetric_mesh(p);
  if (msg != "") {
    throw RuntimeError("Validation of volumetric mesh failed with following error:\n" + msg);
  }
}


} // namespace geometry_utils

} // namespace API
} // namespace MCell