Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80670640
zone.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
Mon, Sep 2, 00:18
Size
13 KB
Mime Type
text/x-c++
Expires
Wed, Sep 4, 00:18 (2 d)
Engine
blob
Format
Raw Data
Handle
20415924
Attached To
rAKA akantu
zone.hh
View Options
/**
* @file zone.hh
*
* @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
*
* @date creation: Fri Jan 04 2013
* @date last modification: Sun Oct 19 2014
*
* @brief contact zone classes
*
* @section LICENSE
*
* Copyright (©) 2014, 2015 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_CZONE_HH__
#define __AKANTU_CZONE_HH__
#include <unordered_set>
#include "aka_common.hh"
#include "surface.hh"
__BEGIN_AKANTU__
//! This class represents a contact element surface
template
<
int
d
>
class
CZone
{
public
:
typedef
CSurface
<
d
>
contact_surface
;
typedef
typename
std
::
list
<
const
contact_surface
*>
surface_list
;
typedef
typename
surface_list
::
const_iterator
surface_iterator
;
typedef
typename
CSurface
<
d
>::
point_type
point_type
;
typedef
typename
CSurface
<
d
>::
bbox_type
bbox_type
;
// typedef std::set<UInt> node_set;
typedef
std
::
unordered_set
<
UInt
>
node_set
;
typedef
typename
node_set
::
iterator
node_iterator
;
typedef
typename
node_set
::
const_iterator
const_node_iterator
;
typedef
std
::
map
<
int
,
node_set
>
bucket_container
;
typedef
typename
bucket_container
::
iterator
bucket_iterator
;
typedef
typename
bucket_container
::
const_iterator
const_bucket_iterator
;
typedef
std
::
set
<
const
Element
*>
element_container
;
typedef
typename
element_container
::
const_iterator
element_iterator
;
// typedef std::map<UInt, element_container> node_element_map;
// typedef typename node_element_map::iterator node_element_iterator;
typedef
std
::
map
<
int
,
element_container
>
elementmap
;
typedef
typename
elementmap
::
iterator
elementmap_iterator
;
typedef
typename
elementmap
::
const_iterator
const_elementmap_iterator
;
private
:
SolidMechanicsModel
&
model_
;
bbox_type
bbox_
;
Real
delta_
[
d
];
//!< grid size length
int
size_
[
d
];
//!< grid dimensions, number of buckets in each dimension
int
index_
[
d
];
//!< indexes to convert to 1D equivalent array
bucket_container
buckets_
;
elementmap
elements_
;
// node_element_map node_element_map_;
surface_list
surfaces_
;
public
:
CZone
()
{}
// clear memory
void
clear
()
{
buckets_
.
clear
();
// node_element_map_.clear();
}
typename
bucket_container
::
size_type
buckets_size
()
const
{
return
buckets_
.
size
();
}
bucket_iterator
buckets_begin
()
{
return
buckets_
.
begin
();
}
bucket_iterator
buckets_end
()
{
return
buckets_
.
end
();
}
const_bucket_iterator
buckets_begin
()
const
{
return
buckets_
.
begin
();
}
const_bucket_iterator
buckets_end
()
const
{
return
buckets_
.
end
();
}
void
erase_bucket
(
bucket_iterator
bit
)
{
buckets_
.
erase
(
bit
);
}
CZone
(
SolidMechanicsModel
&
model
,
const
bbox_type
&
bb
,
const
element_container
&
elems
,
const
contact_surface
&
s1
,
const
contact_surface
&
s2
)
:
model_
(
model
),
bbox_
(
bb
),
delta_
(),
size_
(),
index_
()
{
assert
(
!
elems
.
empty
());
surfaces_
.
push_back
(
&
s1
);
surfaces_
.
push_back
(
&
s2
);
// tolerance
Real
epsilon
=
2
*
std
::
numeric_limits
<
Real
>::
epsilon
();
// loop over elements to get bucket increments
typedef
typename
element_container
::
const_iterator
element_iterator
;
for
(
element_iterator
it
=
elems
.
begin
();
it
!=
elems
.
end
();
++
it
)
{
// get side bounding box
bbox_type
bb
=
getBoundingBox
<
d
>
(
**
it
,
model_
);
// process box to modify delta array
const
point_type
&
m
=
bb
.
min
();
const
point_type
&
M
=
bb
.
max
();
for
(
int
i
=
0
;
i
<
d
;
++
i
)
delta_
[
i
]
=
std
::
max
(
delta_
[
i
],
(
M
[
i
]
-
m
[
i
])
/
3.
+
epsilon
);
}
// get the number of buckets in each direction
const
point_type
&
m
=
bbox_
.
min
();
const
point_type
&
M
=
bbox_
.
max
();
for
(
int
i
=
0
;
i
<
d
;
++
i
)
{
assert
(
delta_
[
i
]
!=
0.
);
size_
[
i
]
=
static_cast
<
int
>
((
M
[
i
]
-
m
[
i
])
/
delta_
[
i
])
+
1
;
}
// get indexes to convert multi-dimensional indexes to a 1D array
index_
[
0
]
=
1
;
for
(
int
i
=
1
;
i
<
d
;
++
i
)
index_
[
i
]
=
index_
[
i
-
1
]
*
size_
[
i
-
1
];
// add elements to buckets
add_elements
(
elems
);
std
::
set
<
UInt
>
nodes_inside
;
Mesh
&
mesh
=
model
.
getMesh
();
const
Array
<
Real
>
&
x
=
model
.
getCurrentPosition
();
// loop over elements
for
(
element_iterator
it
=
elems
.
begin
();
it
!=
elems
.
end
();
++
it
)
{
UInt
nb_nodes
=
mesh
.
getNbNodesPerElement
((
*
it
)
->
type
);
const
Array
<
UInt
>
&
conn
=
mesh
.
getConnectivity
((
*
it
)
->
type
);
// loop over nodes
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
UInt
node
=
conn
((
*
it
)
->
element
,
n
);
point_type
node_coord
=
point_type
(
&
x
(
node
));
// if node is within bounding box
if
(
bbox_
&
node_coord
)
{
nodes_inside
.
insert
(
node
);
// compute index into one-dimensional array
int
coord
[
d
];
for
(
int
i
=
0
;
i
<
d
;
++
i
)
coord
[
i
]
=
static_cast
<
int
>
((
node_coord
[
i
]
-
m
[
i
])
/
delta_
[
i
]);
int
idx
=
compute_index
(
coord
);
// add node to bucket
buckets_
[
idx
].
insert
(
node
);
// // add side to map
// node_element_map_[node].insert(*it);
}
// node within bounding box
}
// loop over nodes
}
// loop over elements
// check nodes
std
::
set
<
UInt
>
node_check
;
for
(
bucket_iterator
it
=
buckets_
.
begin
();
it
!=
buckets_
.
end
();
++
it
)
for
(
node_iterator
nit
=
it
->
second
.
begin
();
nit
!=
it
->
second
.
end
();
++
nit
)
node_check
.
insert
(
*
nit
);
cout
<<
"*** INFO *** Node check after creating buckets passed: "
<<
(
node_check
.
size
()
==
nodes_inside
.
size
())
<<
endl
;
assert
(
node_check
.
size
()
==
nodes_inside
.
size
());
#ifdef DEBUG_CONTACT
cout
<<
"*** INFO *** A total of "
<<
nodes_inside
.
size
()
<<
" nodes lie inside the contact zone."
<<
endl
;
cout
<<
"*** INFO *** A total of "
<<
buckets_
.
size
()
<<
" node buckets were created."
<<
endl
;
#endif
}
int
compute_index
(
int
*
array
)
const
{
int
idx
=
0
;
for
(
int
i
=
0
;
i
<
d
;
++
i
)
idx
+=
index_
[
i
]
*
array
[
i
];
return
idx
;
}
void
decompute_index
(
int
idx
,
int
*
array
,
Int2Type
<
2
>
)
const
{
array
[
0
]
=
idx
%
index_
[
1
];
array
[
1
]
=
idx
/
index_
[
1
];
}
void
decompute_index
(
int
idx
,
int
*
array
,
Int2Type
<
3
>
)
const
{
array
[
0
]
=
(
idx
%
(
index_
[
2
]))
%
index_
[
1
];
array
[
1
]
=
(
idx
%
(
index_
[
2
]))
/
index_
[
1
];
array
[
2
]
=
idx
/
index_
[
2
];
}
void
add_elements
(
const
element_container
&
s
)
{
const
point_type
&
bbm
=
bbox_
.
min
();
// loop over elements
for
(
element_iterator
it
=
s
.
begin
();
it
!=
s
.
end
();
++
it
)
{
bbox_type
elbb
=
getBoundingBox
<
d
>
(
**
it
,
model_
);
const
point_type
&
elmin
=
elbb
.
min
();
const
point_type
&
elmax
=
elbb
.
max
();
int
min
[
d
],
max
[
d
];
for
(
int
i
=
0
;
i
<
d
;
++
i
)
{
min
[
i
]
=
std
::
max
(
0
,
static_cast
<
int
>
((
elmin
[
i
]
-
bbm
[
i
])
/
delta_
[
i
]));
max
[
i
]
=
std
::
min
(
size_
[
i
],
static_cast
<
int
>
((
elmax
[
i
]
-
bbm
[
i
])
/
delta_
[
i
])
+
1
);
}
add_element
(
*
it
,
min
,
max
,
Int2Type
<
d
>
());
}
// loop over elements
}
void
add_element
(
const
Element
*
el
,
int
*
min
,
int
*
max
,
Int2Type
<
2
>
)
{
for
(
int
i
=
min
[
0
];
i
<
max
[
0
];
++
i
)
{
for
(
int
j
=
min
[
1
];
j
<
max
[
1
];
++
j
)
{
int
idxarray
[
d
]
=
{
i
,
j
};
int
idx
=
compute_index
(
idxarray
);
elements_
[
idx
].
insert
(
el
);
}
}
}
void
add_element
(
const
Element
*
el
,
int
*
min
,
int
*
max
,
Int2Type
<
3
>
)
{
for
(
int
i
=
min
[
0
];
i
<
max
[
0
];
++
i
)
{
for
(
int
j
=
min
[
1
];
j
<
max
[
1
];
++
j
)
{
for
(
int
k
=
min
[
2
];
k
<
max
[
2
];
++
k
)
{
int
idxarray
[
d
]
=
{
i
,
j
,
k
};
int
idx
=
compute_index
(
idxarray
);
elements_
[
idx
].
insert
(
el
);
}
}
}
}
element_container
contiguous
(
int
bucket
)
{
element_container
elements
;
collect
(
elements
,
bucket
,
Int2Type
<
d
>
());
return
elements
;
}
bool
in_surface
(
UInt
np
,
const
Element
*
sp
)
{
for
(
surface_iterator
it
=
surfaces_
.
begin
();
it
!=
surfaces_
.
end
();
++
it
)
{
const
contact_surface
&
s
=
**
it
;
if
(
s
.
in_surface
(
np
,
sp
))
return
true
;
}
return
false
;
}
bool
in_element
(
UInt
np
,
const
Element
*
sp
)
{
// loop over nodes of side
Mesh
&
mesh
=
model_
.
getMesh
();
UInt
nb_nodes
=
mesh
.
getNbNodesPerElement
(
sp
->
type
);
const
Array
<
UInt
>
&
conn
=
mesh
.
getConnectivity
(
sp
->
type
);
// loop over nodes
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
UInt
node
=
conn
(
sp
->
element
,
n
);
if
(
node
==
np
)
return
true
;
}
return
false
;
}
void
set_union
(
int
bucket
,
element_container
&
elements
)
const
{
// check if bucket exists
const_elementmap_iterator
it
=
elements_
.
find
(
bucket
);
// if bucket is not found, return since there are no elements to add
if
(
it
==
elements_
.
end
())
return
;
// otherwise add elements
for
(
element_iterator
sit
=
it
->
second
.
begin
();
sit
!=
it
->
second
.
end
();
++
sit
)
elements
.
insert
(
*
sit
);
}
template
<
class
container_type
>
void
collect
(
container_type
&
c
,
int
idx
,
Int2Type
<
2
>
)
const
{
const
int
&
m
=
size_
[
0
];
const
int
&
n
=
size_
[
1
];
int
coord
[
d
];
decompute_index
(
idx
,
coord
,
Int2Type
<
2
>
());
for
(
int
i
=
coord
[
0
]
-
1
;
i
<=
coord
[
0
]
+
1
;
++
i
)
if
(
i
>=
0
&&
i
<
m
)
for
(
int
j
=
coord
[
1
]
-
1
;
j
<=
coord
[
1
]
+
1
;
++
j
)
if
(
j
>=
0
&&
j
<
n
)
{
int
cont
[
d
]
=
{
i
,
j
};
set_union
(
compute_index
(
cont
),
c
);
}
}
template
<
class
container_type
>
void
collect
(
container_type
&
c
,
int
idx
,
Int2Type
<
3
>
)
const
{
const
int
&
m
=
size_
[
0
];
const
int
&
n
=
size_
[
1
];
const
int
&
p
=
size_
[
2
];
int
coord
[
d
];
decompute_index
(
idx
,
coord
,
Int2Type
<
3
>
());
for
(
int
i
=
coord
[
0
]
-
1
;
i
<=
coord
[
0
]
+
1
;
++
i
)
if
(
i
>=
0
&&
i
<
m
)
for
(
int
j
=
coord
[
1
]
-
1
;
j
<=
coord
[
1
]
+
1
;
++
j
)
if
(
j
>=
0
&&
j
<
n
)
for
(
int
k
=
coord
[
2
]
-
1
;
k
<=
coord
[
2
]
+
1
;
++
k
)
if
(
k
>=
0
&&
k
<
p
)
{
int
cont
[
d
]
=
{
i
,
j
,
k
};
set_union
(
compute_index
(
cont
),
c
);
}
}
//! Enable std output
friend
std
::
ostream
&
operator
<<
(
std
::
ostream
&
os
,
const
CZone
&
cz
)
{
Mesh
&
mesh
=
cz
.
model_
.
getMesh
();
const
Array
<
Real
>
&
x
=
cz
.
model_
.
getCurrentPosition
();
int
coord
[
d
];
os
<<
" Contact zone: "
<<
endl
;
os
<<
" origin: "
<<
cz
.
bbox_
.
min
()
<<
endl
;
os
<<
" bounding box: "
<<
cz
.
bbox_
<<
endl
;
os
<<
" delta: ("
;
for
(
int
i
=
0
;
i
<
d
-
1
;
++
i
)
os
<<
cz
.
delta_
[
i
]
<<
","
;
os
<<
cz
.
delta_
[
d
-
1
]
<<
")"
<<
endl
;
os
<<
" size: ("
;
for
(
int
i
=
0
;
i
<
d
-
1
;
++
i
)
os
<<
cz
.
size_
[
i
]
<<
","
;
os
<<
cz
.
size_
[
d
-
1
]
<<
")"
<<
endl
;
os
<<
" number of node buckets: "
<<
cz
.
buckets_
.
size
()
<<
endl
;
std
::
set
<
UInt
>
nodes
;
for
(
const_bucket_iterator
it
=
cz
.
buckets_
.
begin
();
it
!=
cz
.
buckets_
.
end
();
++
it
)
{
os
<<
" node bucket "
<<
it
->
first
;
cz
.
decompute_index
(
it
->
first
,
coord
,
Int2Type
<
d
>
());
os
<<
" ("
<<
coord
[
0
];
for
(
int
i
=
1
;
i
<
d
;
++
i
)
os
<<
","
<<
coord
[
i
];
os
<<
") contains "
<<
it
->
second
.
size
()
<<
" nodes:"
;
for
(
const_node_iterator
nit
=
it
->
second
.
begin
();
nit
!=
it
->
second
.
end
();
++
nit
)
{
os
<<
" "
<<*
nit
<<
" ("
<<
point_type
(
&
x
(
*
nit
))
<<
")"
;
nodes
.
insert
(
*
nit
);
}
os
<<
endl
;
}
os
<<
" there are a total of "
<<
nodes
.
size
()
<<
" nodes in the node buckets."
<<
endl
;
os
<<
" number of element buckets: "
<<
cz
.
elements_
.
size
()
<<
endl
;
for
(
const_elementmap_iterator
it
=
cz
.
elements_
.
begin
();
it
!=
cz
.
elements_
.
end
();
++
it
)
{
os
<<
" element bucket "
<<
it
->
first
;
cz
.
decompute_index
(
it
->
first
,
coord
,
Int2Type
<
d
>
());
os
<<
" ("
<<
coord
[
0
];
for
(
int
i
=
1
;
i
<
d
;
++
i
)
os
<<
","
<<
coord
[
i
];
os
<<
") contains "
<<
it
->
second
.
size
()
<<
" elements:"
<<
endl
;
for
(
element_iterator
sit
=
it
->
second
.
begin
();
sit
!=
it
->
second
.
end
();
++
sit
)
{
os
<<
"
\t\t
"
<<**
sit
;
UInt
nb_nodes
=
mesh
.
getNbNodesPerElement
((
*
sit
)
->
type
);
const
Array
<
UInt
>
&
conn
=
mesh
.
getConnectivity
((
*
sit
)
->
type
);
os
<<
", nodes:"
;
// loop over nodes
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
os
<<
" "
<<
conn
((
*
sit
)
->
element
,
n
)
<<
" ("
<<
point_type
(
&
x
(
conn
((
*
sit
)
->
element
,
n
)))
<<
")"
;
os
<<
endl
;
}
}
return
os
;
}
};
__END_AKANTU__
#endif
/* __AKANTU_CZONE_HH__ */
Event Timeline
Log In to Comment