Homec4science
Diffusion µSpectre (master) Browse

muSpectre/src/common/intersection_octree.cc
eb55b1bae2d8master

Last Change
Raw File

intersection_octree.cc
View Options

/**
* @file intersection_octree.cc
*
* @author Ali Falsafi <ali.falsafi@epfl.ch>
*
* @date May 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.
*/
#include <vector>
#include <array>
#include "common/ccoord_operations.hh"
#include "common/common.hh"
#include "cell/cell_base.hh"
#include "materials/material_base.hh"
#include "common/intersection_octree.hh"
#include "common/intersection_volume_calculator.hh"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template<Dim_t DimS>
Node<DimS>::Node(const Rcoord &new_origin, const Rcoord &new_lengths,
int depth, RootNode_t& root )
:root_node(root),
origin(new_origin),
lengths(new_lengths),
depth(depth),
is_pixel(depth==root.max_depth),
children_no (((is_pixel)? 0 : pow(2, DimS)))
{
this->check_node();
}
template<Dim_t DimS>
Node<DimS>::Node(const CellBase<DimS, DimS>& cell, RootNode_t& root)
:root_node(root),
origin(Rcoord{}),
depth(0),
is_pixel(false),
children_no (((is_pixel)? 0 : pow(2, DimS)))
{
auto cell_resolution = cell.get_domain_resolutions();
auto max_resolution = cell_resolution.at(std::distance(std::max_element(cell_resolution.begin(),
cell_resolution.end())
, cell_resolution.begin()));
auto max_depth = std::ceil(std::log2(max_resolution));
this->lengths = {CcoordOps::get_cube<DimS>(pow(2, max_depth))};
this->check_node();
}
template<Dim_t DimS>
void Node<DimS>::check_node() {
Real total_volume = 1.0, intersected_volume = 0.0, intersection_ratio = 0.0;
for (int i = 0; i < DimS; i++){
total_volume *= this->lengths[i];
}
// this volume should be calculated by CGAL as the intersection volume of the preticipitate and the Node
intersected_volume = 1.0 ;
intersection_volume_calculator<DimS>(this->root_node.preticipitate_vertices,
this->origin, this->lengths);
intersection_ratio = intersected_volume/ total_volume;
if (intersection_ratio == 1.0) {
// all the pixels in this node should be assigned to the precipitate material
Real pix_num = pow(2, (this->root_node.max_depth - this->depth));
Ccoord origin_point, pixels_number;
for (int i = 0; i < DimS; i++ ) {
origin_point[i] = this->root_node.origin.at(i) / this->root_node.pixel_lengths.at(i);
pixels_number[i] = pix_num;
}
CcoordOps::Pixels<DimS> pixels (pixels_number, origin_point);
for (auto && pix:pixels){
if (this->root_node.cell.is_inside(origin_point)){
this->root_node.intersected_pixels.push_back(pix);
this->root_node.intersection_ratios.push_back(intersection_ratio);
}
}
}else if (intersection_ratio > 0) {
this->split_node(intersection_ratio);
}
}
template<Dim_t DimS>
void Node<DimS>::split_node(Real intersection_ratio){
if (this->is_pixel) {
Ccoord pixel;
for (int i = 0 ; i < DimS ; i++){
pixel[i] = this->origin[i] / this->root_node.pixel_lengths[i];
}
this->root_node.intersected_pixels.push_back(pixel);
this->root_node.intersection_ratios.push_back(intersection_ratio);
} else {
this->divide_node();
}
}
template<Dim_t DimS>
void Node<DimS>::divide_node(){
Rcoord new_origin, new_length;
this->children.reserve(children_no);
CcoordOps::Pixels<DimS> sub_nodes (CcoordOps::get_cube<DimS>(2));
for (auto && sub_node: sub_nodes ) {
for (int i = 0 ; i< DimS ; i++){
new_length[i] = this->lengths[i] *0.5;
new_origin[i] = this->origin[i] + sub_node[i] * new_length[i];
}
this->children.emplace_back(new_origin, new_length,this->depth+1, this->root_node);
}
}
template<Dim_t DimS>
RootNode<DimS>::RootNode(CellBase<DimS, DimS>& cell,
std::vector<Rcoord> vert_precipitate)
//Calling parent constructin method simply by argumnets of (coordinates of origin, 2^depth_max in each dirextion
//,0 ,*this)
:Node<DimS>(Rcoord{},
CcoordOps::get_cube<DimS>(pow(2,
std::ceil(std::log2(this->cell.get_domain_resolutions().
at(std::distance(std::max_element
(cell.get_domain_resolutions().begin(),
cell.get_domain_resolutions().end())
,cell.get_domain_resolutions().begin())))))),
0, *this),
cell(cell),
cell_length (cell.get_domain_lengths()),
pixel_lengths (cell.get_pixel_lengths()),
cell_resolution (cell.get_domain_resolutions()),
max_resolution (cell_resolution.at(std::distance(std::max_element(cell_resolution.begin(),
cell_resolution.end())
,cell_resolution.begin()))),
max_depth (std::ceil(std::log2(this->max_resolution))),
preticipitate_vertices {vert_precipitate}
{
this->check_node();
}
template<Dim_t DimS>
std::vector<Ccoord_t<DimS>> RootNode<DimS>::get_intersected_pixels(){
return this->intersected_pixels;
}
template<Dim_t DimS>
std::vector<Real> RootNode<DimS>::get_intersection_ratios(){
return this->intersection_ratios;
}
template class RootNode<threeD>;
template class RootNode<twoD>;
template class Node<threeD>;
template class Node<twoD>;
} //muSpctre
c4science · Help