Homec4science
Diffusion µSpectre (feat/repair_intersection_octree) Browse

muSpectre/src/common/intersection_octree.cc
80eb1f984de2feat/repair_intersectio…

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 "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"
#include <vector>
#include <array>
#include <algorithm>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template<Dim_t DimS>
Node<DimS>::Node(const Rcoord &new_origin, const Ccoord &new_lengths,
int depth, RootNode_t& root, bool is_root )
:root_node(root),
origin(new_origin),
Clengths(new_lengths),
depth(depth),
is_pixel((depth == root.max_depth)),
children_no (((is_pixel)? 0 : pow(2, DimS)))
{
for (int i = 0 ; i < DimS ; i++){
this->Rlengths[i] = this->Clengths[i] * this->root_node.pixel_lengths[i];
}
if (not is_root){
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;
constexpr Real machine_prescision = 1e-12;
for (int i = 0; i < DimS; i++){
total_volume *= this->Rlengths[i];
}
// this volume should be calculated by CGAL as the intersection volume of the precipitate and the Node
intersected_volume = total_volume ;
intersected_volume =
intersection_volume_calculator<DimS>(this->root_node.precipitate_vertices,
this->origin, this->Rlengths);
intersection_ratio = intersected_volume / total_volume;
if (intersection_ratio > (1.0 - machine_prescision)) {
// 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] = std::round(this->origin[i] / this->root_node.pixel_lengths[i]);
pixels_number[i] = pix_num;
}
CcoordOps::Pixels<DimS> pixels(pixels_number, origin_point);
for (auto && pix:pixels){
this->root_node.intersected_pixels.push_back(pix);
this->root_node.intersection_ratios.push_back(1.0);
}
}else if (intersection_ratio > machine_prescision) {
this->split_node(intersection_ratio);
}
}
/* ---------------------------------------------------------------------- */
template<Dim_t DimS>
void Node<DimS>::split_node(Real intersection_ratio){
if (this->depth == this->root_node.max_depth) {
Ccoord pixel;
for (int i = 0 ; i < DimS ; i++){
pixel[i] = std::round(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;
Ccoord 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] = std::round(this->Clengths[i] *0.5);
new_origin[i] = this->origin[i] + sub_node[i] * Rlengths[i] * 0.5;
}
this->children.emplace_back(new_origin, new_length, this->depth+1, this->root_node, false);
}
}
/* ---------------------------------------------------------------------- */
template<Dim_t DimS>
RootNode<DimS>::RootNode(CellBase<DimS, DimS>& cell,
std::vector<Rcoord> vert_precipitate)
/*Calling parent constructing method simply by argumnets of (coordinates of origin, 2^depth_max in each direction, 0 ,*this)*/
:Node<DimS>(Rcoord{},
CcoordOps::get_cube<DimS, Dim_t>(pow(2,
std::ceil(std::log2(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, true),
cell(cell),
cell_length (cell.get_domain_lengths()),
pixel_lengths (cell.get_pixel_lengths()),
cell_resolution (cell.get_domain_resolutions()),
max_resolution (this->cell_resolution.at(std::distance(std::max_element(this->cell_resolution.begin(),
this->cell_resolution.end())
,this->cell_resolution.begin()))),
max_depth (std::ceil(std::log2(this->max_resolution))),
precipitate_vertices (vert_precipitate)
{
for (int i = 0 ; i < DimS ; i++){
this->Rlengths[i] = this->Clengths[i] * this->root_node.pixel_lengths[i];
}
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<Dim_t DimS>
void RootNode<DimS>:: check_root_node()
{
for (int i = 0 ; i < DimS ; i++){
this->Rlengths[i] = this->Clengths[i] * this->root_node.pixel_lengths[i];
}
this->check_node();
}
/* ---------------------------------------------------------------------- */
template class RootNode<threeD>;
template class RootNode<twoD>;
template class Node<threeD>;
template class Node<twoD>;
} //muSpctre
c4science · Help