Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80410917
mesh_segment_intersector_tmpl.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Aug 31, 14:21
Size
10 KB
Mime Type
text/x-c++
Expires
Mon, Sep 2, 14:21 (2 d)
Engine
blob
Format
Raw Data
Handle
20003761
Attached To
rAKA akantu
mesh_segment_intersector_tmpl.hh
View Options
/**
* @file mesh_segment_intersector_tmpl.hh
*
* @author Lucas Frerot <lucas.frerot@epfl.ch>
* @author Clement Roux <clement.roux@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Wed Apr 29 2015
* @date last modification: Tue Feb 20 2018
*
* @brief Computation of mesh intersection with segments
*
* @section LICENSE
*
* Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MESH_SEGMENT_INTERSECTOR_TMPL_HH__
#define __AKANTU_MESH_SEGMENT_INTERSECTOR_TMPL_HH__
#include "aka_common.hh"
#include "mesh_geom_common.hh"
#include "tree_type_helper.hh"
namespace
akantu
{
template
<
UInt
dim
,
ElementType
type
>
MeshSegmentIntersector
<
dim
,
type
>::
MeshSegmentIntersector
(
Mesh
&
mesh
,
Mesh
&
result_mesh
)
:
parent_type
(
mesh
),
result_mesh
(
result_mesh
),
current_physical_name
()
{
this
->
intersection_points
=
new
Array
<
Real
>
(
0
,
dim
);
this
->
constructData
();
}
template
<
UInt
dim
,
ElementType
type
>
MeshSegmentIntersector
<
dim
,
type
>::~
MeshSegmentIntersector
()
{}
template
<
UInt
dim
,
ElementType
type
>
void
MeshSegmentIntersector
<
dim
,
type
>::
computeIntersectionQuery
(
const
K
::
Segment_3
&
query
)
{
AKANTU_DEBUG_IN
();
result_mesh
.
addConnectivityType
(
_segment_2
,
_not_ghost
);
result_mesh
.
addConnectivityType
(
_segment_2
,
_ghost
);
std
::
list
<
result_type
>
result_list
;
std
::
set
<
std
::
pair
<
K
::
Segment_3
,
UInt
>
,
segmentPairsLess
>
segment_set
;
this
->
factory
.
getTree
().
all_intersections
(
query
,
std
::
back_inserter
(
result_list
));
this
->
computeSegments
(
result_list
,
segment_set
,
query
);
// Arrays for storing nodes and connectivity
Array
<
Real
>
&
nodes
=
result_mesh
.
getNodes
();
Array
<
UInt
>
&
connectivity
=
result_mesh
.
getConnectivity
(
_segment_2
);
// Arrays for storing associated element and physical name
bool
valid_elemental_data
=
true
;
Array
<
Element
>
*
associated_element
=
NULL
;
Array
<
std
::
string
>
*
associated_physical_name
=
NULL
;
try
{
associated_element
=
&
result_mesh
.
getData
<
Element
>
(
"associated_element"
,
_segment_2
);
associated_physical_name
=
&
result_mesh
.
getData
<
std
::
string
>
(
"physical_names"
,
_segment_2
);
}
catch
(
debug
::
Exception
&
e
)
{
valid_elemental_data
=
false
;
}
std
::
set
<
pair_type
,
segmentPairsLess
>::
iterator
it
=
segment_set
.
begin
(),
end
=
segment_set
.
end
();
// Loop over the segment pairs
for
(;
it
!=
end
;
++
it
)
{
if
(
!
it
->
first
.
is_degenerate
())
{
Vector
<
UInt
>
segment_connectivity
(
2
);
segment_connectivity
(
0
)
=
result_mesh
.
getNbNodes
();
segment_connectivity
(
1
)
=
result_mesh
.
getNbNodes
()
+
1
;
connectivity
.
push_back
(
segment_connectivity
);
// Copy nodes
Vector
<
Real
>
source
(
dim
),
target
(
dim
);
for
(
UInt
j
=
0
;
j
<
dim
;
j
++
)
{
source
(
j
)
=
it
->
first
.
source
()[
j
];
target
(
j
)
=
it
->
first
.
target
()[
j
];
}
nodes
.
push_back
(
source
);
nodes
.
push_back
(
target
);
// Copy associated element info
if
(
valid_elemental_data
)
{
associated_element
->
push_back
(
Element
{
type
,
it
->
second
,
_not_ghost
});
associated_physical_name
->
push_back
(
current_physical_name
);
}
}
}
AKANTU_DEBUG_OUT
();
}
template
<
UInt
dim
,
ElementType
type
>
void
MeshSegmentIntersector
<
dim
,
type
>::
computeMeshQueryIntersectionPoint
(
const
K
::
Segment_3
&
/*query*/
,
UInt
/*nb_old_nodes*/
)
{
AKANTU_ERROR
(
"The method: computeMeshQueryIntersectionPoint has not "
"been implemented in class MeshSegmentIntersector!"
);
}
template
<
UInt
dim
,
ElementType
type
>
void
MeshSegmentIntersector
<
dim
,
type
>::
buildResultFromQueryList
(
const
std
::
list
<
K
::
Segment_3
>
&
query_list
)
{
AKANTU_DEBUG_IN
();
this
->
computeIntersectionQueryList
(
query_list
);
AKANTU_DEBUG_OUT
();
}
template
<
UInt
dim
,
ElementType
type
>
void
MeshSegmentIntersector
<
dim
,
type
>::
computeSegments
(
const
std
::
list
<
result_type
>
&
intersections
,
std
::
set
<
pair_type
,
segmentPairsLess
>
&
segments
,
const
K
::
Segment_3
&
query
)
{
AKANTU_DEBUG_IN
();
/*
* Number of intersections = 0 means
*
* - query is completely outside mesh
* - query is completely inside primitive
*
* We try to determine the case and still construct the segment list
*/
if
(
intersections
.
size
()
==
0
)
{
// We look at all the primitives intersected by two rays
// If there is one primitive in common, then query is inside
// that primitive
K
::
Ray_3
ray1
(
query
.
source
(),
query
.
target
());
K
::
Ray_3
ray2
(
query
.
target
(),
query
.
source
());
std
::
set
<
UInt
>
ray1_results
,
ray2_results
;
this
->
factory
.
getTree
().
all_intersected_primitives
(
ray1
,
std
::
inserter
(
ray1_results
,
ray1_results
.
begin
()));
this
->
factory
.
getTree
().
all_intersected_primitives
(
ray2
,
std
::
inserter
(
ray2_results
,
ray2_results
.
begin
()));
bool
inside_primitive
=
false
;
UInt
primitive_id
=
0
;
std
::
set
<
UInt
>::
iterator
ray2_it
=
ray2_results
.
begin
(),
ray2_end
=
ray2_results
.
end
();
// Test if first list contains an element of second list
for
(;
ray2_it
!=
ray2_end
&&
!
inside_primitive
;
++
ray2_it
)
{
if
(
ray1_results
.
find
(
*
ray2_it
)
!=
ray1_results
.
end
())
{
inside_primitive
=
true
;
primitive_id
=
*
ray2_it
;
}
}
if
(
inside_primitive
)
{
segments
.
insert
(
std
::
make_pair
(
query
,
primitive_id
));
}
}
else
{
typename
std
::
list
<
result_type
>::
const_iterator
it
=
intersections
.
begin
(),
end
=
intersections
.
end
();
for
(;
it
!=
end
;
++
it
)
{
UInt
el
=
(
*
it
)
->
second
;
// Result of intersection is a segment
if
(
const
K
::
Segment_3
*
segment
=
boost
::
get
<
K
::
Segment_3
>
(
&
((
*
it
)
->
first
)))
{
// Check if the segment was alread created
segments
.
insert
(
std
::
make_pair
(
*
segment
,
el
));
}
// Result of intersection is a point
else
if
(
const
K
::
Point_3
*
point
=
boost
::
get
<
K
::
Point_3
>
(
&
((
*
it
)
->
first
)))
{
// We only want to treat points differently if we're in 3D with Tetra4
// elements This should be optimized by compilator
if
(
dim
==
3
&&
type
==
_tetrahedron_4
)
{
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
TreeTypeHelper
<
Triangle
<
K
>
,
K
>::
container_type
facets
;
const
Array
<
Real
>
&
nodes
=
this
->
mesh
.
getNodes
();
Array
<
UInt
>::
const_vector_iterator
connectivity_vec
=
this
->
mesh
.
getConnectivity
(
type
).
begin
(
nb_nodes_per_element
);
const
Vector
<
UInt
>
&
el_connectivity
=
connectivity_vec
[
el
];
Matrix
<
Real
>
node_coordinates
(
dim
,
nb_nodes_per_element
);
for
(
UInt
i
=
0
;
i
<
nb_nodes_per_element
;
i
++
)
for
(
UInt
j
=
0
;
j
<
dim
;
j
++
)
node_coordinates
(
j
,
i
)
=
nodes
(
el_connectivity
(
i
),
j
);
this
->
factory
.
addPrimitive
(
node_coordinates
,
el
,
facets
);
// Local tree
TreeTypeHelper
<
Triangle
<
K
>
,
K
>::
tree
*
local_tree
=
new
TreeTypeHelper
<
Triangle
<
K
>
,
K
>::
tree
(
facets
.
begin
(),
facets
.
end
());
// Compute local intersections (with current element)
std
::
list
<
result_type
>
local_intersections
;
local_tree
->
all_intersections
(
query
,
std
::
back_inserter
(
local_intersections
));
bool
out_point_found
=
false
;
typename
std
::
list
<
result_type
>::
const_iterator
local_it
=
local_intersections
.
begin
(),
local_end
=
local_intersections
.
end
();
for
(;
local_it
!=
local_end
;
++
local_it
)
{
if
(
const
K
::
Point_3
*
local_point
=
boost
::
get
<
K
::
Point_3
>
(
&
((
*
local_it
)
->
first
)))
{
if
(
!
comparePoints
(
*
point
,
*
local_point
))
{
K
::
Segment_3
seg
(
*
point
,
*
local_point
);
segments
.
insert
(
std
::
make_pair
(
seg
,
el
));
out_point_found
=
true
;
}
}
}
if
(
!
out_point_found
)
{
TreeTypeHelper
<
Triangle
<
K
>
,
K
>::
point_type
a
(
node_coordinates
(
0
,
0
),
node_coordinates
(
1
,
0
),
node_coordinates
(
2
,
0
)),
b
(
node_coordinates
(
0
,
1
),
node_coordinates
(
1
,
1
),
node_coordinates
(
2
,
1
)),
c
(
node_coordinates
(
0
,
2
),
node_coordinates
(
1
,
2
),
node_coordinates
(
2
,
2
)),
d
(
node_coordinates
(
0
,
3
),
node_coordinates
(
1
,
3
),
node_coordinates
(
2
,
3
));
K
::
Tetrahedron_3
tetra
(
a
,
b
,
c
,
d
);
const
K
::
Point_3
*
inside_point
=
NULL
;
if
(
tetra
.
has_on_bounded_side
(
query
.
source
())
&&
!
tetra
.
has_on_boundary
(
query
.
source
()))
inside_point
=
&
query
.
source
();
else
if
(
tetra
.
has_on_bounded_side
(
query
.
target
())
&&
!
tetra
.
has_on_boundary
(
query
.
target
()))
inside_point
=
&
query
.
target
();
if
(
inside_point
)
{
K
::
Segment_3
seg
(
*
inside_point
,
*
point
);
segments
.
insert
(
std
::
make_pair
(
seg
,
el
));
}
}
delete
local_tree
;
}
}
}
}
AKANTU_DEBUG_OUT
();
}
}
// namespace akantu
#endif
// __AKANTU_MESH_SEGMENT_INTERSECTOR_TMPL_HH__
Event Timeline
Log In to Comment