Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78005753
test_cohesive_parallel_buildfragments.cc
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 17, 19:28
Size
13 KB
Mime Type
text/x-c
Expires
Mon, Aug 19, 19:28 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
19852440
Attached To
rAKA akantu
test_cohesive_parallel_buildfragments.cc
View Options
/**
* @file test_cohesive_parallel_buildfragments.cc
*
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
*
* @brief Test to build fragments in parallel
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
*/
/* -------------------------------------------------------------------------- */
#include <limits>
#include <fstream>
#include <iostream>
#include <algorithm>
#include <functional>
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model_cohesive.hh"
#include "material_cohesive.hh"
#include "fragment_manager.hh"
/* -------------------------------------------------------------------------- */
using
namespace
akantu
;
void
verticalInsertionLimit
(
SolidMechanicsModelCohesive
&
);
void
displaceElements
(
SolidMechanicsModelCohesive
&
,
const
Real
,
const
Real
);
bool
isInertiaEqual
(
const
Vector
<
Real
>
&
,
const
Vector
<
Real
>
&
);
void
rotateArray
(
Array
<
Real
>
&
array
,
Real
angle
);
UInt
getNbBigFragments
(
FragmentManager
&
,
UInt
);
const
UInt
spatial_dimension
=
3
;
const
UInt
total_nb_fragment
=
4
;
const
Real
rotation_angle
=
M_PI
/
4.
;
const
Real
global_tolerance
=
1.e-9
;
int
main
(
int
argc
,
char
*
argv
[])
{
initialize
(
"material.dat"
,
argc
,
argv
);
Math
::
setTolerance
(
global_tolerance
);
Mesh
mesh
(
spatial_dimension
);
StaticCommunicator
&
comm
=
StaticCommunicator
::
getStaticCommunicator
();
Int
psize
=
comm
.
getNbProc
();
Int
prank
=
comm
.
whoAmI
();
akantu
::
MeshPartition
*
partition
=
NULL
;
if
(
prank
==
0
)
{
// Read the mesh
mesh
.
read
(
"mesh.msh"
);
/// partition the mesh
MeshUtils
::
purifyMesh
(
mesh
);
partition
=
new
MeshPartitionScotch
(
mesh
,
spatial_dimension
);
partition
->
partitionate
(
psize
);
}
SolidMechanicsModelCohesive
model
(
mesh
);
model
.
initParallel
(
partition
,
NULL
,
true
);
delete
partition
;
/// model initialization
model
.
initFull
(
SolidMechanicsModelCohesiveOptions
(
_explicit_lumped_mass
,
true
));
mesh
.
computeBoundingBox
();
Real
L
=
mesh
.
getUpperBounds
()(
0
)
-
mesh
.
getLowerBounds
()(
0
);
Real
h
=
mesh
.
getUpperBounds
()(
1
)
-
mesh
.
getLowerBounds
()(
1
);
Real
rho
=
model
.
getMaterial
(
"bulk"
).
getParam
<
Real
>
(
"rho"
);
Real
theoretical_mass
=
L
*
h
*
h
*
rho
;
Real
frag_theo_mass
=
theoretical_mass
/
total_nb_fragment
;
UInt
nb_element
=
mesh
.
getNbElement
(
spatial_dimension
,
_not_ghost
,
_ek_regular
);
comm
.
allReduce
(
&
nb_element
,
1
,
_so_sum
);
UInt
nb_element_per_fragment
=
nb_element
/
total_nb_fragment
;
FragmentManager
fragment_manager
(
model
);
fragment_manager
.
computeAllData
();
getNbBigFragments
(
fragment_manager
,
nb_element_per_fragment
+
1
);
model
.
setBaseName
(
"extrinsic"
);
model
.
addDumpFieldVector
(
"displacement"
);
model
.
addDumpField
(
"velocity"
);
model
.
addDumpField
(
"stress"
);
model
.
addDumpField
(
"partitions"
);
model
.
addDumpField
(
"fragments"
);
model
.
addDumpField
(
"fragments mass"
);
model
.
addDumpField
(
"moments of inertia"
);
model
.
addDumpField
(
"elements per fragment"
);
model
.
dump
();
model
.
setBaseNameToDumper
(
"cohesive elements"
,
"cohesive_elements_test"
);
model
.
addDumpFieldVectorToDumper
(
"cohesive elements"
,
"displacement"
);
model
.
addDumpFieldToDumper
(
"cohesive elements"
,
"damage"
);
model
.
dump
(
"cohesive elements"
);
/// set check facets
verticalInsertionLimit
(
model
);
model
.
assembleMassLumped
();
model
.
synchronizeBoundaries
();
/// impose initial displacement
Array
<
Real
>
&
displacement
=
model
.
getDisplacement
();
Array
<
Real
>
&
velocity
=
model
.
getVelocity
();
const
Array
<
Real
>
&
position
=
mesh
.
getNodes
();
UInt
nb_nodes
=
mesh
.
getNbNodes
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
displacement
(
n
,
0
)
=
position
(
n
,
0
)
*
0.1
;
velocity
(
n
,
0
)
=
position
(
n
,
0
);
}
rotateArray
(
mesh
.
getNodes
(),
rotation_angle
);
// rotateArray(displacement, rotation_angle);
// rotateArray(velocity, rotation_angle);
model
.
updateResidual
();
model
.
checkCohesiveStress
();
model
.
dump
();
model
.
dump
(
"cohesive elements"
);
const
Array
<
Real
>
&
fragment_mass
=
fragment_manager
.
getMass
();
const
Array
<
Real
>
&
fragment_center
=
fragment_manager
.
getCenterOfMass
();
Real
el_size
=
L
/
total_nb_fragment
;
Real
lim
=
-
L
/
2
+
el_size
*
0.99
;
/// define theoretical inertia moments
Vector
<
Real
>
small_frag_inertia
(
spatial_dimension
);
small_frag_inertia
(
0
)
=
frag_theo_mass
*
(
h
*
h
+
h
*
h
)
/
12.
;
small_frag_inertia
(
1
)
=
frag_theo_mass
*
(
el_size
*
el_size
+
h
*
h
)
/
12.
;
small_frag_inertia
(
2
)
=
frag_theo_mass
*
(
el_size
*
el_size
+
h
*
h
)
/
12.
;
std
::
sort
(
small_frag_inertia
.
storage
(),
small_frag_inertia
.
storage
()
+
spatial_dimension
,
std
::
greater
<
Real
>
());
const
Array
<
Real
>
&
inertia_moments
=
fragment_manager
.
getMomentsOfInertia
();
Array
<
Real
>::
const_iterator
<
Vector
<
Real
>
>
inertia_moments_begin
=
inertia_moments
.
begin
(
spatial_dimension
);
/// displace one fragment each time
for
(
UInt
frag
=
1
;
frag
<=
total_nb_fragment
;
++
frag
)
{
if
(
prank
==
0
)
std
::
cout
<<
"Generating fragment: "
<<
frag
<<
std
::
endl
;
fragment_manager
.
computeAllData
();
/// check number of big fragments
UInt
nb_big_fragment
=
getNbBigFragments
(
fragment_manager
,
nb_element_per_fragment
+
1
);
model
.
dump
();
model
.
dump
(
"cohesive elements"
);
if
(
frag
<
total_nb_fragment
)
{
if
(
nb_big_fragment
!=
1
)
{
AKANTU_DEBUG_ERROR
(
"The number of big fragments is wrong: "
<<
nb_big_fragment
);
return
EXIT_FAILURE
;
}
}
else
{
if
(
nb_big_fragment
!=
0
)
{
AKANTU_DEBUG_ERROR
(
"The number of big fragments is wrong: "
<<
nb_big_fragment
);
return
EXIT_FAILURE
;
}
}
/// check number of fragments
UInt
nb_fragment_num
=
fragment_manager
.
getNbFragment
();
if
(
nb_fragment_num
!=
frag
)
{
AKANTU_DEBUG_ERROR
(
"The number of fragments is wrong! Numerical: "
<<
nb_fragment_num
<<
" Theoretical: "
<<
frag
);
return
EXIT_FAILURE
;
}
/// check mass computation
if
(
frag
<
total_nb_fragment
)
{
Real
total_mass
=
0.
;
UInt
small_fragments
=
0
;
for
(
UInt
f
=
0
;
f
<
nb_fragment_num
;
++
f
)
{
const
Vector
<
Real
>
&
current_inertia
=
inertia_moments_begin
[
f
];
if
(
Math
::
are_float_equal
(
fragment_mass
(
f
,
0
),
frag_theo_mass
))
{
/// check center of mass
if
(
fragment_center
(
f
,
0
)
>
(
L
*
frag
/
total_nb_fragment
-
L
/
2
))
{
AKANTU_DEBUG_ERROR
(
"Fragment center is wrong!"
);
return
EXIT_FAILURE
;
}
/// check moment of inertia
if
(
!
isInertiaEqual
(
current_inertia
,
small_frag_inertia
))
{
AKANTU_DEBUG_ERROR
(
"Inertia moments are wrong"
);
return
EXIT_FAILURE
;
}
++
small_fragments
;
total_mass
+=
frag_theo_mass
;
}
else
{
/// check the moment of inertia for the biggest fragment
Real
big_frag_mass
=
frag_theo_mass
*
(
total_nb_fragment
-
frag
+
1
);
Real
big_frag_size
=
el_size
*
(
total_nb_fragment
-
frag
+
1
);
Vector
<
Real
>
big_frag_inertia
(
spatial_dimension
);
big_frag_inertia
(
0
)
=
big_frag_mass
*
(
h
*
h
+
h
*
h
)
/
12.
;
big_frag_inertia
(
1
)
=
big_frag_mass
*
(
big_frag_size
*
big_frag_size
+
h
*
h
)
/
12.
;
big_frag_inertia
(
2
)
=
big_frag_mass
*
(
big_frag_size
*
big_frag_size
+
h
*
h
)
/
12.
;
std
::
sort
(
big_frag_inertia
.
storage
(),
big_frag_inertia
.
storage
()
+
spatial_dimension
,
std
::
greater
<
Real
>
());
if
(
!
isInertiaEqual
(
current_inertia
,
big_frag_inertia
))
{
AKANTU_DEBUG_ERROR
(
"Inertia moments are wrong"
);
return
EXIT_FAILURE
;
}
}
}
if
(
small_fragments
!=
nb_fragment_num
-
1
)
{
AKANTU_DEBUG_ERROR
(
"The number of small fragments is wrong!"
);
return
EXIT_FAILURE
;
}
if
(
!
Math
::
are_float_equal
(
total_mass
,
small_fragments
*
frag_theo_mass
))
{
AKANTU_DEBUG_ERROR
(
"The mass of small fragments is wrong!"
);
return
EXIT_FAILURE
;
}
}
/// displace fragments
rotateArray
(
mesh
.
getNodes
(),
-
rotation_angle
);
// rotateArray(displacement, -rotation_angle);
displaceElements
(
model
,
lim
,
el_size
*
2
);
rotateArray
(
mesh
.
getNodes
(),
rotation_angle
);
// rotateArray(displacement, rotation_angle);
model
.
updateResidual
();
lim
+=
el_size
;
}
model
.
dump
();
model
.
dump
(
"cohesive elements"
);
/// check centers
const
Array
<
Real
>
&
fragment_velocity
=
fragment_manager
.
getVelocity
();
Real
initial_position
=
-
L
/
2.
+
el_size
/
2.
;
for
(
UInt
frag
=
0
;
frag
<
total_nb_fragment
;
++
frag
)
{
Real
theoretical_center
=
initial_position
+
el_size
*
frag
;
UInt
f_index
=
0
;
while
(
f_index
<
total_nb_fragment
&&
!
Math
::
are_float_equal
(
fragment_center
(
f_index
,
0
),
theoretical_center
))
++
f_index
;
if
(
f_index
==
total_nb_fragment
)
{
AKANTU_DEBUG_ERROR
(
"The fragments' center is wrong!"
);
return
EXIT_FAILURE
;
}
f_index
=
0
;
while
(
f_index
<
total_nb_fragment
&&
!
Math
::
are_float_equal
(
fragment_velocity
(
f_index
,
0
),
theoretical_center
))
++
f_index
;
if
(
f_index
==
total_nb_fragment
)
{
AKANTU_DEBUG_ERROR
(
"The fragments' velocity is wrong!"
);
return
EXIT_FAILURE
;
}
}
finalize
();
if
(
prank
==
0
)
std
::
cout
<<
"OK: test_cohesive_buildfragments was passed!"
<<
std
::
endl
;
return
EXIT_SUCCESS
;
}
void
verticalInsertionLimit
(
SolidMechanicsModelCohesive
&
model
)
{
UInt
spatial_dimension
=
model
.
getSpatialDimension
();
const
Mesh
&
mesh_facets
=
model
.
getMeshFacets
();
const
Array
<
Real
>
&
position
=
mesh_facets
.
getNodes
();
for
(
ghost_type_t
::
iterator
gt
=
ghost_type_t
::
begin
();
gt
!=
ghost_type_t
::
end
();
++
gt
)
{
GhostType
ghost_type
=
*
gt
;
Mesh
::
type_iterator
it
=
mesh_facets
.
firstType
(
spatial_dimension
-
1
,
ghost_type
);
Mesh
::
type_iterator
end
=
mesh_facets
.
lastType
(
spatial_dimension
-
1
,
ghost_type
);
for
(;
it
!=
end
;
++
it
)
{
ElementType
type
=
*
it
;
Array
<
bool
>
&
check_facets
=
model
.
getElementInserter
().
getCheckFacets
(
type
,
ghost_type
);
const
Array
<
UInt
>
&
connectivity
=
mesh_facets
.
getConnectivity
(
type
,
ghost_type
);
UInt
nb_nodes_per_facet
=
connectivity
.
getNbComponent
();
for
(
UInt
f
=
0
;
f
<
check_facets
.
getSize
();
++
f
)
{
if
(
!
check_facets
(
f
))
continue
;
UInt
nb_aligned_nodes
=
1
;
Real
first_node_pos
=
position
(
connectivity
(
f
,
0
),
0
);
for
(;
nb_aligned_nodes
<
nb_nodes_per_facet
;
++
nb_aligned_nodes
)
{
Real
other_node_pos
=
position
(
connectivity
(
f
,
nb_aligned_nodes
),
0
);
if
(
!
Math
::
are_float_equal
(
first_node_pos
,
other_node_pos
))
break
;
}
if
(
nb_aligned_nodes
!=
nb_nodes_per_facet
)
{
check_facets
(
f
)
=
false
;
}
}
}
}
}
void
displaceElements
(
SolidMechanicsModelCohesive
&
model
,
const
Real
lim
,
const
Real
amount
)
{
UInt
spatial_dimension
=
model
.
getSpatialDimension
();
Array
<
Real
>
&
displacement
=
model
.
getDisplacement
();
Mesh
&
mesh
=
model
.
getMesh
();
UInt
nb_nodes
=
mesh
.
getNbNodes
();
Array
<
bool
>
displaced
(
nb_nodes
);
displaced
.
clear
();
Vector
<
Real
>
barycenter
(
spatial_dimension
);
for
(
ghost_type_t
::
iterator
gt
=
ghost_type_t
::
begin
();
gt
!=
ghost_type_t
::
end
();
++
gt
)
{
GhostType
ghost_type
=
*
gt
;
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
,
ghost_type
);
Mesh
::
type_iterator
end
=
mesh
.
lastType
(
spatial_dimension
,
ghost_type
);
for
(;
it
!=
end
;
++
it
)
{
ElementType
type
=
*
it
;
const
Array
<
UInt
>
&
connectivity
=
mesh
.
getConnectivity
(
type
,
ghost_type
);
UInt
nb_element
=
connectivity
.
getSize
();
UInt
nb_nodes_per_element
=
connectivity
.
getNbComponent
();
Array
<
UInt
>::
const_vector_iterator
conn_el
=
connectivity
.
begin
(
nb_nodes_per_element
);
for
(
UInt
el
=
0
;
el
<
nb_element
;
++
el
)
{
mesh
.
getBarycenter
(
el
,
type
,
barycenter
.
storage
(),
ghost_type
);
if
(
barycenter
(
0
)
<
lim
)
{
const
Vector
<
UInt
>
&
conn
=
conn_el
[
el
];
for
(
UInt
n
=
0
;
n
<
nb_nodes_per_element
;
++
n
)
{
UInt
node
=
conn
(
n
);
if
(
!
displaced
(
node
))
{
displacement
(
node
,
0
)
-=
amount
;
displaced
(
node
)
=
true
;
}
}
}
}
}
}
}
bool
isInertiaEqual
(
const
Vector
<
Real
>
&
a
,
const
Vector
<
Real
>
&
b
)
{
UInt
nb_terms
=
a
.
size
();
UInt
equal_terms
=
0
;
while
(
equal_terms
<
nb_terms
&&
std
::
abs
(
a
(
equal_terms
)
-
b
(
equal_terms
))
/
a
(
equal_terms
)
<
Math
::
getTolerance
())
++
equal_terms
;
return
equal_terms
==
nb_terms
;
}
void
rotateArray
(
Array
<
Real
>
&
array
,
Real
angle
)
{
UInt
spatial_dimension
=
array
.
getNbComponent
();
Real
rotation_values
[]
=
{
std
::
cos
(
angle
),
std
::
sin
(
angle
),
0
,
-
std
::
sin
(
angle
),
std
::
cos
(
angle
),
0
,
0
,
0
,
1
};
Matrix
<
Real
>
rotation
(
rotation_values
,
spatial_dimension
,
spatial_dimension
);
RVector
displaced_node
(
spatial_dimension
);
Array
<
Real
>::
vector_iterator
node_it
=
array
.
begin
(
spatial_dimension
);
Array
<
Real
>::
vector_iterator
node_end
=
array
.
end
(
spatial_dimension
);
for
(;
node_it
!=
node_end
;
++
node_it
)
{
displaced_node
.
mul
<
false
>
(
rotation
,
*
node_it
);
*
node_it
=
displaced_node
;
}
}
UInt
getNbBigFragments
(
FragmentManager
&
fragment_manager
,
UInt
minimum_nb_elements
)
{
fragment_manager
.
computeNbElementsPerFragment
();
const
Array
<
UInt
>
&
nb_elements_per_fragment
=
fragment_manager
.
getNbElementsPerFragment
();
UInt
nb_fragment
=
fragment_manager
.
getNbFragment
();
UInt
nb_big_fragment
=
0
;
for
(
UInt
frag
=
0
;
frag
<
nb_fragment
;
++
frag
)
{
if
(
nb_elements_per_fragment
(
frag
)
>=
minimum_nb_elements
)
{
++
nb_big_fragment
;
}
}
return
nb_big_fragment
;
}
Event Timeline
Log In to Comment