33 : surface_region_(surface_region) {
43 if (surface_vec(2)==0) {
44 THROW( ExcSurfaceProjection() << EI_Message(
"Vertical plane of surface is forbidden!") );
50 ta /= arma::norm(ta, 2);
52 tb /= arma::norm(tb, 2);
66 if (region_set.size() == 0)
67 THROW( RegionDB::ExcUnknownSet() << RegionDB::EI_Label(surface_region) );
68 for (
auto reg : region_set) {
70 THROW( ExcSurfaceProjection()
71 << EI_Message(
"Dimension of surface region " + surface_region +
" must be 2!") );
79 if (ele.region().is_in_region_set(region_set)) {
83 project_node(0) = projection(0); project_node(1) = projection(1);
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);
93 a_mat.col(0) = *ele.node(1) - *ele.node(0);
94 a_mat.col(1) = *ele.node(2) - *ele.node(0);
98 b_vecs_.push_back( *ele.node(0) );
103 if ( boxes.size() == 0) {
104 THROW( ExcSurfaceProjection()
105 << EI_Message(
"Region " + surface_region +
" contains no boundary element! Probably bulk region was set.") );
114 double distance = std::numeric_limits<double>::max();
115 bool found_surface_projection =
false;
119 project_point.subvec(0,1) =
m_ * point;
120 project_point(2) = 0;
126 if ( (x(0)>=0) && (x(1)>=0) && (x(0)+x(1)<=1) ) {
127 double new_distance = -x(2);
128 if ( (new_distance>=0) && (new_distance<distance) ) distance = new_distance;
129 found_surface_projection =
true;
133 if (!found_surface_projection) {
134 double snap_dist = std::numeric_limits<double>::max();
154 auto proj_3d = M.submat(0,0,2,1) * proj_ref +
b_vecs_[*
it];
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;
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_) );
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
static const unsigned int default_leaf_size_limit
Default leaf size limit.
void add_boxes(const std::vector< BoundingBox > &boxes)
void find_point(const Space< 3 >::Point &point, std::vector< unsigned int > &result_list, bool full_list=false) const
Bounding box in 3d ambient space.
void expand(const Point &point)
const RegionDB & region_db() const
Range< ElementAccessor< 3 > > elements_range() const
Returns range of mesh elements.
BCMesh * bc_mesh() const override
Implement MeshBase::bc_mesh(), getter of boundary mesh.
double global_snap_radius() const
Maximal distance of observe point from Mesh relative to its size.
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
static BaryPoint local_to_bary(const LocalPoint &lp)
Converts from local to barycentric coordinates.
static BaryPoint clip(const BaryPoint &barycentric)
RegionSet get_region_set(const std::string &set_name) const
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.
ArmaMat< double, N, M > mat
Implementation of range helper class.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.