Homec4science
Diffusion µSpectre (intersection_octree) Browse

muSpectre/src/common/intersection_volume_calculator.hh
3b9c93d008c2intersection_octree

Last Change
Raw File

intersection_volume_calculator.hh
View Options

/**
* @file intersection_volume_calculator.hh
*
* @author Ali Falsafi <ali.falsafi@epfl.ch>
*
* @date 04 June 2018
*
* @brief common operations on pixel addressing
*
* Copyright © 2018 Ali Falsafi
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre 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
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef INTERSECTION_VOLUME_CALCULATOR_H
#define INTERSECTION_VOLUME_CALCULATOR_H
#include <CGAL/Polyhedron_incremental_builder_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/Nef_polyhedron_3.h>
#include <CGAL/IO/Nef_polyhedron_iostream_3.h>
#include <CGAL/Polyhedron_items_with_id_3.h>
#include <CGAL/Aff_transformation_3.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/Surface_mesh.h>
#include "common/ccoord_operations.hh"
#include "common/common.hh"
#include <vector>
#include <fstream>
#include <math.h>
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef CGAL::Polyhedron_3<Kernel, CGAL::Polyhedron_items_with_id_3> Polyhedron;
typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron;
typedef Polyhedron::HalfedgeDS HalfedgeDS;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Vector_3 Vector_3;
typedef Kernel::Tetrahedron_3 Tetrahedron;
typedef Polyhedron::Vertex_iterator Vertex_iterator;
typedef Polyhedron::Facet_iterator Facet_iterator;
typedef Polyhedron::Halfedge_around_facet_circulator Hafc;
typedef CGAL::Aff_transformation_3<Kernel> Transformation_mat;
using Dim_t = int;
using Real = double;
Polyhedron node_polyhedron_generator(std::array<Real, 3> origin, std::array<Real, 3> lengths) {
std::vector<Point_3> node_vertices ;
node_vertices.push_back(Point_3(origin[0] , origin[1] , origin[2] ));
node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] , origin[2] ));
node_vertices.push_back(Point_3(origin[0] , origin[1] + lengths[1], origin[2] ));
node_vertices.push_back(Point_3(origin[0] , origin[1] , origin[2] + lengths[2]));
node_vertices.push_back(Point_3(- origin[0] , origin[1] + lengths[1], origin[2] + lengths[2]));
node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] , origin[2] + lengths[2]));
node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] + lengths[1], origin[2] ));
node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] + lengths[1], origin[2] + lengths[2]));
Polyhedron node_poly;
CGAL::convex_hull_3(node_vertices.begin(), node_vertices.end(), node_poly);
return node_poly;
}
Polyhedron convex_polyhedron_generator(std::vector<std::array<Real, 3>> convex_poly_vertices) {
std::vector<Point_3> poly_vertices ;
for (auto && point:convex_poly_vertices ){
poly_vertices.push_back(Point_3(point[0], point[1], point[2]));
}
Polyhedron convex_poly;
CGAL::convex_hull_3(poly_vertices.begin(), poly_vertices.end(), convex_poly);
return convex_poly;
}
template <Dim_t DimS, class T, class U>
class Correction{
public:
U correct_origin (T array);
U correct_length (T array);
U correct_vector (T vector);
};
template<>
class Correction <2, std::array<Real,2>, std::array<Real,3>> {
public:
std::array<Real,3> correct_origin (std::array<Real,2> array){
return std::array<Real,3>{array.at(0),array.at(1),0.0};
}
std::array<Real,3> correct_length (std::array<Real,2> array){
return std::array<Real,3>{array.at(0),array.at(1),1.0};
}
};
template<>
class Correction <3, std::array<Real,3>, std::array<Real,3>> {
public:
std::array<Real,3> correct_origin (std::array<Real,3> array){
return array;
}
std::array<Real,3> correct_length (std::array<Real,3> array){
return array;
}
};
template<>
class Correction <2, std::vector<std::array<Real,2>>, std::vector<std::array<Real,3>>> {
public:
std::vector<std::array<Real, 3>> correct_vector(std::vector<std::array<Real, 2>> vertices){
std::vector<std::array<Real, 3>> corrected_convex_poly_vertices;
for (auto && vertice : vertices){
corrected_convex_poly_vertices.push_back({vertice.at(0), vertice.at(1), 0.0});
corrected_convex_poly_vertices.push_back({vertice.at(0), vertice.at(1), 1.0});
}
return corrected_convex_poly_vertices;
}
};
template<>
class Correction <3, std::vector<std::array<Real,3>>, std::vector<std::array<Real,3>>> {
public:
std::vector<std::array<Real, 3>> correct_vector(std::vector<std::array<Real, 3>> vertices){
std::vector<std::array<Real, 3>> corrected_convex_poly_vertices;
return vertices;
}
};
template <Dim_t DimS>
Real intersection_volume_calculator(std::vector<std::array<Real, DimS>> convex_poly_vertices,
std::array<Real, DimS> origin, std::array<Real, DimS> lengths){
Correction<DimS, std::array<Real,DimS>, std::array<Real,3>> array_correction;
Correction<DimS, std::vector<std::array<Real,DimS>>, std::vector<std::array<Real,3>>> vector_correction;
std::vector<std::array<Real,3>> corrected_convex_poly_vertices (vector_correction.correct_vector(convex_poly_vertices));
std::array<Real,3> corrected_origin(array_correction.correct_origin(origin));
std::array<Real,3> corrected_lengths(array_correction.correct_length(lengths));
std::vector<Point_3> result_vertex ;
std::vector<std::size_t> a_result_facet;
Vector_3 result_vector ;
Point_3 center_point = Point_3( 0.0, 0.0, 0.0) ;
Point_3 origin_coordinate_point (0,0,0), tet_vertices[4];
Real return_volume = 0.0;
Polyhedron node_poly, convex_poly, result_poly;
Tetrahedron tetra;
node_poly = node_polyhedron_generator(corrected_origin, corrected_lengths);
convex_poly = convex_polyhedron_generator(corrected_convex_poly_vertices);
Nef_polyhedron convex_poly_nef, node_poly_nef, result_nef ;
convex_poly_nef = Nef_polyhedron (convex_poly);
node_poly_nef = Nef_polyhedron (node_poly);
//computing the intersection:
result_nef = convex_poly_nef * node_poly_nef ;
if(result_nef.is_simple() and not result_nef.is_empty()) {
result_nef.convert_to_Polyhedron(result_poly);
//extracting vertex
int v_id = 0 ;
for ( Vertex_iterator v = result_poly.vertices_begin(); v != result_poly.vertices_end(); ++v){
result_vertex.push_back(v->point());
result_vector = v->point() - origin_coordinate_point;
center_point = center_point + result_vector;
v->id() = v_id++;
}
//calculating the center point of the resultant polyhedron:
Transformation_mat scale (1,0,0,0,0,1,0,0,0,0,1,0,v_id);
center_point = scale(center_point);
tet_vertices[0] = center_point;
int f_id = 0;
for ( Facet_iterator f = result_poly.facets_begin(); f != result_poly.facets_end(); ++f){
a_result_facet.clear();
f->id() = f_id++;
Hafc circ = f->facet_begin();
//assignig corresponding vertex to facets:
do {
a_result_facet.push_back(circ-> vertex()-> id());
} while ( ++circ != f->facet_begin());
//making tetreahedrals to calculate volumes:
for(int ii = 0 ; ii<3; ii++) tet_vertices[ii+1] = result_vertex[a_result_facet[ii]];
//calculating volumes:
tetra = Tetrahedron (tet_vertices[0], tet_vertices[1], tet_vertices[2], tet_vertices[3]);
return_volume += CGAL::to_double(tetra.volume());
}
}
return return_volume;
}
#endif //INTERSECTION_VOLUME_CALCULATOR
c4science · Help