Flow123d  master-f44eb46
surface_depth.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file surface_depth.cc
15  * @brief
16  */
17 
18 #include <vector>
19 #include <armadillo>
20 #include <limits>
21 
22 #include "fields/surface_depth.hh"
23 #include "mesh/accessors.hh"
24 #include "mesh/mesh.h"
25 #include "mesh/bc_mesh.hh"
26 #include "mesh/ref_element.hh"
27 #include "system/sys_profiler.hh"
28 #include "mesh/accessors.hh"
29 #include "mesh/range_wrapper.hh"
30 
31 
32 SurfaceDepth::SurfaceDepth(const Mesh *mesh, std::string surface_region, std::string surface_direction)
33 : surface_region_(surface_region) {
34  this->create_projection_matrix( arma::vec3(surface_direction) );
35  this->construct_bih_tree( const_cast<Mesh*>(mesh), surface_region);
38 }
39 
40 
42 {
43  if (surface_vec(2)==0) {
44  THROW( ExcSurfaceProjection() << EI_Message("Vertical plane of surface is forbidden!") );
45  }
46  surface_norm_vec_ = surface_vec / arma::norm(surface_vec, 2); // normalize to size == 1
47 
48  arma::vec3 ex("1 0 0"), ey("0 1 0");
49  arma::vec3 ta = ex - arma::dot(ex, surface_norm_vec_) * surface_norm_vec_;
50  ta /= arma::norm(ta, 2); // normalize
51  arma::vec3 tb = ey - arma::dot(ey, surface_norm_vec_) * surface_norm_vec_;
52  tb /= arma::norm(tb, 2); // normalize
53  m_ = arma::mat(2,3);
54  m_.row(0) = ta.t();
55  m_.row(1) = tb.t();
56 }
57 
58 
59 void SurfaceDepth::construct_bih_tree(Mesh *mesh, std::string surface_region)
60 {
62  inv_projection_.clear();
63  b_vecs_.clear();
64 
65  RegionSet region_set = mesh->region_db().get_region_set(surface_region);
66  if (region_set.size() == 0)
67  THROW( RegionDB::ExcUnknownSet() << RegionDB::EI_Label(surface_region) );
68  for (auto reg : region_set) {
69  if (reg.dim() != 2)
70  THROW( ExcSurfaceProjection()
71  << EI_Message("Dimension of surface region " + surface_region + " must be 2!") );
72  }
73 
74  // make element boxes
75  unsigned int i=0;
76  unsigned int i_node;
77  arma::vec3 project_node("0 0 0");
78  for( auto ele : mesh->bc_mesh()->elements_range() ) {
79  if (ele.region().is_in_region_set(region_set)) {
80  ASSERT_EQ(ele->n_nodes(), 3);
81 
82  arma::vec projection = m_ * (*ele.node(0));
83  project_node(0) = projection(0); project_node(1) = projection(1);
84  BoundingBox bb(project_node);
85  for(i_node=1; i_node<ele->n_nodes(); i_node++) {
86  arma::vec project_coords = m_ * (*ele.node(i_node));
87  project_node(0) = project_coords(0); project_node(1) = project_coords(1);
88  bb.expand(project_node);
89  }
90  boxes.push_back(bb);
91 
92  arma::mat a_mat(3,3);
93  a_mat.col(0) = *ele.node(1) - *ele.node(0);
94  a_mat.col(1) = *ele.node(2) - *ele.node(0);
95  a_mat.col(2) = surface_norm_vec_;
96 
97  inv_projection_.push_back( a_mat.i() );
98  b_vecs_.push_back( *ele.node(0) );
99  }
100  i++;
101  }
102 
103  if ( boxes.size() == 0) {
104  THROW( ExcSurfaceProjection()
105  << EI_Message("Region " + surface_region + " contains no boundary element! Probably bulk region was set.") );
106  }
107 
108  bih_tree_.add_boxes( boxes );
110 }
111 
113 {
114  double distance = std::numeric_limits<double>::max();
115  bool found_surface_projection = false;
116  arma::vec3 project_point;
117  arma::vec3 x;
118 
119  project_point.subvec(0,1) = m_ * point;
120  project_point(2) = 0;
121 
122  searched_elements_.clear();
123  bih_tree_.find_point(project_point, searched_elements_);
125  prepare_distance_solve( (*it), point, x );
126  if ( (x(0)>=0) && (x(1)>=0) && (x(0)+x(1)<=1) ) { // point is in triangle
127  double new_distance = -x(2);
128  if ( (new_distance>=0) && (new_distance<distance) ) distance = new_distance;
129  found_surface_projection = true;
130  }
131  }
132 
133  if (!found_surface_projection) {
134  double snap_dist = std::numeric_limits<double>::max();
135  arma::vec3 proj_to_surface_plane;
136 
138  // check snap distance point - triangle
139 
140  // compute `x` - local element coordinates of `point`
141  prepare_distance_solve( (*it), point, x );
142 
143  // project `point` to the surface plane of the found element
144  proj_to_surface_plane = point - x(2)*surface_norm_vec_;
145 
146  // get local x,y coordinates inside the found element
147  arma::vec local_point = x.subvec(0,1);
148  // transform to barycentric coords and clip
149  auto bary_point = RefElement<2>::local_to_bary(local_point);
150  auto clip_point = RefElement<2>::clip(bary_point);
151  auto proj_ref = RefElement<2>::bary_to_local(clip_point);
152  // get the element map and map from unit to real (like in mapping)
153  arma::mat M = inv_projection_[*it].i();
154  auto proj_3d = M.submat(0,0,2,1) * proj_ref + b_vecs_[*it];
155 
156  // difference between snap and projected point
157  double new_snap_dist = arma::norm(proj_3d - proj_to_surface_plane, 2);
158  if (new_snap_dist<snap_dist) {
159  snap_dist = new_snap_dist;
160  distance = -x(2);
161  }
162  }
163 
164  if (snap_dist > projection_search_radius_) {
165  THROW(ExcTooLargeSnapDistance() << EI_Xcoord(proj_to_surface_plane(0)) << EI_Ycoord(proj_to_surface_plane(1))
166  << EI_Zcoord(proj_to_surface_plane(2)) << EI_SnapDistance(snap_dist) << EI_RegionName(surface_region_) );
167  }
168  }
169 
170  return distance;
171 }
172 
173 void SurfaceDepth::prepare_distance_solve(unsigned int elem_idx, arma::vec3 &point, arma::vec3 &x)
174 {
175  x = inv_projection_[elem_idx] * (point - b_vecs_[elem_idx]);
176 }
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
void construct()
Definition: bih_tree.cc:55
static const unsigned int default_leaf_size_limit
Default leaf size limit.
Definition: bih_tree.hh:45
void add_boxes(const std::vector< BoundingBox > &boxes)
Definition: bih_tree.cc:43
void find_point(const Space< 3 >::Point &point, std::vector< unsigned int > &result_list, bool full_list=false) const
Definition: bih_tree.cc:287
Bounding box in 3d ambient space.
Definition: bounding_box.hh:53
void expand(const Point &point)
const RegionDB & region_db() const
Definition: mesh.h:175
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
Definition: mesh.cc:1174
Definition: mesh.h:362
BCMesh * bc_mesh() const override
Implement MeshBase::bc_mesh(), getter of boundary mesh.
Definition: mesh.h:567
double global_snap_radius() const
Maximal distance of observe point from Mesh relative to its size.
Definition: mesh.cc:1077
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
Definition: ref_element.cc:201
static BaryPoint local_to_bary(const LocalPoint &lp)
Converts from local to barycentric coordinates.
Definition: ref_element.cc:187
static BaryPoint clip(const BaryPoint &barycentric)
Definition: ref_element.cc:332
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:328
void create_projection_matrix(arma::vec3 surface_vec)
Create projection matrix m_.
void construct_bih_tree(Mesh *mesh, std::string surface_region)
Construct BIH tree above surface region of given name.
double projection_search_radius_
Projection search radius given from global_snap_radius of mesh.
std::vector< arma::vec3 > b_vecs_
vector of b-vectors (right side of equation system) of elements above which BIH tree is created
double compute_distance(arma::vec3 point)
Compute distance of point from given surface region (was set in constructor)
std::vector< unsigned int > searched_elements_
Surface region name (need for exception)
std::string surface_region_
Surface region name (need for exception)
arma::mat m_
projection matrix
void prepare_distance_solve(unsigned int elem_idx, arma::vec3 &point, arma::vec3 &x)
Precompute solve of distance in transformed coordinations of surface element plane.
arma::vec3 surface_norm_vec_
normal vector of surface plane
std::vector< arma::mat > inv_projection_
vector of inversion projection matrices of elements above which BIH tree is created
SurfaceDepth(const Mesh *mesh, std::string surface_region, std::string surface_direction)
BIHTree bih_tree_
Tree of mesh elements.
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
ArmaMat< double, N, M > mat
Definition: armor.hh:888
ArmaVec< double, N > vec
Definition: armor.hh:885
Implementation of range helper class.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.