Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78256390
solid_mechanics_model.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
Mon, Aug 19, 11:02
Size
47 KB
Mime Type
text/x-c
Expires
Wed, Aug 21, 11:02 (2 d)
Engine
blob
Format
Raw Data
Handle
19955920
Attached To
rAKA akantu
solid_mechanics_model.cc
View Options
/**
* @file solid_mechanics_model.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Fri Sep 19 2014
*
* @brief Implementation of the SolidMechanicsModel class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_math.hh"
#include "aka_common.hh"
#include "solid_mechanics_model.hh"
#include "group_manager_inline_impl.cc"
#include "dumpable_inline_impl.hh"
#include "integration_scheme_2nd_order.hh"
#include "element_group.hh"
#include "static_communicator.hh"
#include "dof_synchronizer.hh"
#include "element_group.hh"
#include <cmath>
#ifdef AKANTU_USE_IOHELPER
# include "dumper_field.hh"
# include "dumper_paraview.hh"
# include "dumper_homogenizing_field.hh"
# include "dumper_material_internal_field.hh"
# include "dumper_elemental_field.hh"
# include "dumper_material_padders.hh"
# include "dumper_element_partition.hh"
# include "dumper_iohelper.hh"
#endif
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
const
SolidMechanicsModelOptions
default_solid_mechanics_model_options
(
_explicit_lumped_mass
,
false
);
/* -------------------------------------------------------------------------- */
/**
* A solid mechanics model need a mesh and a dimension to be created. the model
* by it self can not do a lot, the good init functions should be called in
* order to configure the model depending on what we want to do.
*
* @param mesh mesh representing the model we want to simulate
* @param dim spatial dimension of the problem, if dim = 0 (default value) the
* dimension of the problem is assumed to be the on of the mesh
* @param id an id to identify the model
*/
SolidMechanicsModel
::
SolidMechanicsModel
(
Mesh
&
mesh
,
UInt
dim
,
const
ID
&
id
,
const
MemoryID
&
memory_id
)
:
Model
(
mesh
,
dim
,
id
,
memory_id
),
BoundaryCondition
<
SolidMechanicsModel
>
(),
time_step
(
NAN
),
f_m2a
(
1.0
),
mass_matrix
(
NULL
),
velocity_damping_matrix
(
NULL
),
stiffness_matrix
(
NULL
),
jacobian_matrix
(
NULL
),
material_index
(
"material index"
,
id
),
material_local_numbering
(
"material local numbering"
,
id
),
material_selector
(
new
DefaultMaterialSelector
(
material_index
)),
is_default_material_selector
(
true
),
integrator
(
NULL
),
increment_flag
(
false
),
synch_parallel
(
NULL
),
are_materials_instantiated
(
false
)
{
AKANTU_DEBUG_IN
();
createSynchronizerRegistry
(
this
);
registerFEEngineObject
<
MyFEEngineType
>
(
"SolidMechanicsFEEngine"
,
mesh
,
spatial_dimension
);
this
->
displacement
=
NULL
;
this
->
mass
=
NULL
;
this
->
velocity
=
NULL
;
this
->
acceleration
=
NULL
;
this
->
force
=
NULL
;
this
->
internal_force
=
NULL
;
this
->
residual
=
NULL
;
this
->
blocked_dofs
=
NULL
;
this
->
increment
=
NULL
;
this
->
increment_acceleration
=
NULL
;
this
->
dof_synchronizer
=
NULL
;
this
->
previous_displacement
=
NULL
;
materials
.
clear
();
mesh
.
registerEventHandler
(
*
this
);
#if defined(AKANTU_USE_IOHELPER)
this
->
mesh
.
registerDumper
<
DumperParaview
>
(
"paraview_all"
,
id
,
true
);
this
->
mesh
.
addDumpMesh
(
mesh
,
spatial_dimension
,
_not_ghost
,
_ek_regular
);
#endif
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel
::~
SolidMechanicsModel
()
{
AKANTU_DEBUG_IN
();
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
delete
*
mat_it
;
}
materials
.
clear
();
delete
integrator
;
delete
synch_parallel
;
if
(
is_default_material_selector
)
{
delete
material_selector
;
material_selector
=
NULL
;
}
AKANTU_DEBUG_OUT
();
}
void
SolidMechanicsModel
::
setTimeStep
(
Real
time_step
)
{
this
->
time_step
=
time_step
;
#if defined(AKANTU_USE_IOHELPER)
this
->
mesh
.
getDumper
().
setTimeStep
(
time_step
);
#endif
}
/* -------------------------------------------------------------------------- */
/* Initialization */
/* -------------------------------------------------------------------------- */
/**
* This function groups many of the initialization in on function. For most of
* basics case the function should be enough. The functions initialize the
* model, the internal vectors, set them to 0, and depending on the parameters
* it also initialize the explicit or implicit solver.
*
* @param material_file the file containing the materials to use
* @param method the analysis method wanted. See the akantu::AnalysisMethod for
* the different possibilities
*/
void
SolidMechanicsModel
::
initFull
(
const
ModelOptions
&
options
)
{
Model
::
initFull
(
options
);
const
SolidMechanicsModelOptions
&
smm_options
=
dynamic_cast
<
const
SolidMechanicsModelOptions
&>
(
options
);
method
=
smm_options
.
analysis_method
;
// initialize the vectors
initArrays
();
// set the initial condition to 0
force
->
clear
();
velocity
->
clear
();
acceleration
->
clear
();
displacement
->
clear
();
// initialize pbc
if
(
pbc_pair
.
size
()
!=
0
)
initPBC
();
// initialize the time integration schemes
switch
(
method
)
{
case
_explicit_lumped_mass:
initExplicit
();
break
;
case
_explicit_consistent_mass:
initSolver
();
initExplicit
();
break
;
case
_implicit_dynamic:
initImplicit
(
true
);
break
;
case
_static:
initImplicit
(
false
);
break
;
default
:
AKANTU_EXCEPTION
(
"analysis method not recognised by SolidMechanicsModel"
);
break
;
}
// initialize the materials
if
(
this
->
parser
->
getLastParsedFile
()
!=
""
)
{
instantiateMaterials
();
}
if
(
!
smm_options
.
no_init_materials
)
{
initMaterials
();
}
if
(
increment_flag
)
initBC
(
*
this
,
*
displacement
,
*
increment
,
*
force
);
else
initBC
(
*
this
,
*
displacement
,
*
force
);
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
initParallel
(
MeshPartition
*
partition
,
DataAccessor
*
data_accessor
)
{
AKANTU_DEBUG_IN
();
if
(
data_accessor
==
NULL
)
data_accessor
=
this
;
synch_parallel
=
&
createParallelSynch
(
partition
,
data_accessor
);
synch_registry
->
registerSynchronizer
(
*
synch_parallel
,
_gst_material_id
);
synch_registry
->
registerSynchronizer
(
*
synch_parallel
,
_gst_smm_mass
);
synch_registry
->
registerSynchronizer
(
*
synch_parallel
,
_gst_smm_stress
);
synch_registry
->
registerSynchronizer
(
*
synch_parallel
,
_gst_smm_boundary
);
synch_registry
->
registerSynchronizer
(
*
synch_parallel
,
_gst_for_dump
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
initFEEngineBoundary
()
{
FEEngine
&
fem_boundary
=
getFEEngineBoundary
();
fem_boundary
.
initShapeFunctions
(
_not_ghost
);
fem_boundary
.
initShapeFunctions
(
_ghost
);
fem_boundary
.
computeNormalsOnControlPoints
(
_not_ghost
);
fem_boundary
.
computeNormalsOnControlPoints
(
_ghost
);
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
initExplicit
(
AnalysisMethod
analysis_method
)
{
AKANTU_DEBUG_IN
();
//in case of switch from implicit to explicit
if
(
!
this
->
isExplicit
())
method
=
analysis_method
;
if
(
integrator
)
delete
integrator
;
integrator
=
new
CentralDifference
();
UInt
nb_nodes
=
acceleration
->
getSize
();
UInt
nb_degree_of_freedom
=
acceleration
->
getNbComponent
();
std
::
stringstream
sstr
;
sstr
<<
id
<<
":increment_acceleration"
;
increment_acceleration
=
&
(
alloc
<
Real
>
(
sstr
.
str
(),
nb_nodes
,
nb_degree_of_freedom
,
Real
()));
AKANTU_DEBUG_OUT
();
}
void
SolidMechanicsModel
::
initArraysPreviousDisplacment
()
{
AKANTU_DEBUG_IN
();
SolidMechanicsModel
::
setIncrementFlagOn
();
UInt
nb_nodes
=
mesh
.
getNbNodes
();
std
::
stringstream
sstr_disp_t
;
sstr_disp_t
<<
id
<<
":previous_displacement"
;
previous_displacement
=
&
(
alloc
<
Real
>
(
sstr_disp_t
.
str
(),
nb_nodes
,
spatial_dimension
,
0.
));
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
/**
* Allocate all the needed vectors. By default their are not necessarily set to
* 0
*
*/
void
SolidMechanicsModel
::
initArrays
()
{
AKANTU_DEBUG_IN
();
UInt
nb_nodes
=
mesh
.
getNbNodes
();
std
::
stringstream
sstr_disp
;
sstr_disp
<<
id
<<
":displacement"
;
// std::stringstream sstr_mass; sstr_mass << id << ":mass";
std
::
stringstream
sstr_velo
;
sstr_velo
<<
id
<<
":velocity"
;
std
::
stringstream
sstr_acce
;
sstr_acce
<<
id
<<
":acceleration"
;
std
::
stringstream
sstr_forc
;
sstr_forc
<<
id
<<
":force"
;
std
::
stringstream
sstr_ifor
;
sstr_forc
<<
id
<<
":internal_force"
;
std
::
stringstream
sstr_resi
;
sstr_resi
<<
id
<<
":residual"
;
std
::
stringstream
sstr_boun
;
sstr_boun
<<
id
<<
":blocked_dofs"
;
displacement
=
&
(
alloc
<
Real
>
(
sstr_disp
.
str
(),
nb_nodes
,
spatial_dimension
,
REAL_INIT_VALUE
));
// mass = &(alloc<Real>(sstr_mass.str(), nb_nodes, spatial_dimension, 0));
velocity
=
&
(
alloc
<
Real
>
(
sstr_velo
.
str
(),
nb_nodes
,
spatial_dimension
,
REAL_INIT_VALUE
));
acceleration
=
&
(
alloc
<
Real
>
(
sstr_acce
.
str
(),
nb_nodes
,
spatial_dimension
,
REAL_INIT_VALUE
));
force
=
&
(
alloc
<
Real
>
(
sstr_forc
.
str
(),
nb_nodes
,
spatial_dimension
,
REAL_INIT_VALUE
));
internal_force
=
&
(
alloc
<
Real
>
(
sstr_ifor
.
str
(),
nb_nodes
,
spatial_dimension
,
REAL_INIT_VALUE
));
residual
=
&
(
alloc
<
Real
>
(
sstr_resi
.
str
(),
nb_nodes
,
spatial_dimension
,
REAL_INIT_VALUE
));
blocked_dofs
=
&
(
alloc
<
bool
>
(
sstr_boun
.
str
(),
nb_nodes
,
spatial_dimension
,
false
));
std
::
stringstream
sstr_curp
;
sstr_curp
<<
id
<<
":current_position"
;
current_position
=
&
(
alloc
<
Real
>
(
sstr_curp
.
str
(),
0
,
spatial_dimension
,
REAL_INIT_VALUE
));
for
(
UInt
g
=
_not_ghost
;
g
<=
_ghost
;
++
g
)
{
GhostType
gt
=
(
GhostType
)
g
;
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
,
gt
,
_ek_not_defined
);
Mesh
::
type_iterator
end
=
mesh
.
lastType
(
spatial_dimension
,
gt
,
_ek_not_defined
);
for
(;
it
!=
end
;
++
it
)
{
UInt
nb_element
=
mesh
.
getNbElement
(
*
it
,
gt
);
material_index
.
alloc
(
nb_element
,
1
,
*
it
,
gt
);
material_local_numbering
.
alloc
(
nb_element
,
1
,
*
it
,
gt
);
}
}
dof_synchronizer
=
new
DOFSynchronizer
(
mesh
,
spatial_dimension
);
dof_synchronizer
->
initLocalDOFEquationNumbers
();
dof_synchronizer
->
initGlobalDOFEquationNumbers
();
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*
*/
void
SolidMechanicsModel
::
initModel
()
{
/// \todo add the current position as a parameter to initShapeFunctions for
/// large deformation
getFEEngine
().
initShapeFunctions
(
_not_ghost
);
getFEEngine
().
initShapeFunctions
(
_ghost
);
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
initPBC
()
{
Model
::
initPBC
();
registerPBCSynchronizer
();
// as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves
std
::
map
<
UInt
,
UInt
>::
iterator
it
=
pbc_pair
.
begin
();
std
::
map
<
UInt
,
UInt
>::
iterator
end
=
pbc_pair
.
end
();
UInt
dim
=
mesh
.
getSpatialDimension
();
while
(
it
!=
end
)
{
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
(
*
blocked_dofs
)((
*
it
).
first
,
i
)
=
true
;
++
it
;
}
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
registerPBCSynchronizer
(){
PBCSynchronizer
*
synch
=
new
PBCSynchronizer
(
pbc_pair
);
synch_registry
->
registerSynchronizer
(
*
synch
,
_gst_smm_uv
);
synch_registry
->
registerSynchronizer
(
*
synch
,
_gst_smm_mass
);
synch_registry
->
registerSynchronizer
(
*
synch
,
_gst_smm_res
);
synch_registry
->
registerSynchronizer
(
*
synch
,
_gst_for_dump
);
changeLocalEquationNumberForPBC
(
pbc_pair
,
mesh
.
getSpatialDimension
());
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
updateCurrentPosition
()
{
AKANTU_DEBUG_IN
();
UInt
nb_nodes
=
mesh
.
getNbNodes
();
current_position
->
resize
(
nb_nodes
);
Real
*
current_position_val
=
current_position
->
storage
();
Real
*
position_val
=
mesh
.
getNodes
().
storage
();
Real
*
displacement_val
=
displacement
->
storage
();
/// compute current_position = initial_position + displacement
memcpy
(
current_position_val
,
position_val
,
nb_nodes
*
spatial_dimension
*
sizeof
(
Real
));
for
(
UInt
n
=
0
;
n
<
nb_nodes
*
spatial_dimension
;
++
n
)
{
*
current_position_val
++
+=
*
displacement_val
++
;
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
initializeUpdateResidualData
()
{
AKANTU_DEBUG_IN
();
UInt
nb_nodes
=
mesh
.
getNbNodes
();
residual
->
resize
(
nb_nodes
);
/// copy the forces in residual for boundary conditions
memcpy
(
residual
->
storage
(),
force
->
storage
(),
nb_nodes
*
spatial_dimension
*
sizeof
(
Real
));
// start synchronization
synch_registry
->
asynchronousSynchronize
(
_gst_smm_uv
);
synch_registry
->
waitEndSynchronize
(
_gst_smm_uv
);
updateCurrentPosition
();
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
/* Explicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* This function compute the second member of the motion equation. That is to
* say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given
* by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$
* F_{int} = \int_{\Omega} N \sigma d\Omega@f$
*
*/
void
SolidMechanicsModel
::
updateResidual
(
bool
need_initialize
)
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_INFO
(
"Assemble the internal forces"
);
// f = f_ext - f_int
// f = f_ext
if
(
need_initialize
)
initializeUpdateResidualData
();
AKANTU_DEBUG_INFO
(
"Compute local stresses"
);
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
mat
.
computeAllStresses
(
_not_ghost
);
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* ------------------------------------------------------------------------ */
/* Computation of the non local part */
synch_registry
->
asynchronousSynchronize
(
_gst_mnl_for_average
);
AKANTU_DEBUG_INFO
(
"Compute non local stresses for local elements"
);
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
mat
.
computeAllNonLocalStresses
(
_not_ghost
);
}
AKANTU_DEBUG_INFO
(
"Wait distant non local stresses"
);
synch_registry
->
waitEndSynchronize
(
_gst_mnl_for_average
);
AKANTU_DEBUG_INFO
(
"Compute non local stresses for ghosts elements"
);
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
mat
.
computeAllNonLocalStresses
(
_ghost
);
}
#endif
/* ------------------------------------------------------------------------ */
/* assembling the forces internal */
// communicate the stress
AKANTU_DEBUG_INFO
(
"Send data for residual assembly"
);
synch_registry
->
asynchronousSynchronize
(
_gst_smm_stress
);
AKANTU_DEBUG_INFO
(
"Assemble residual for local elements"
);
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
mat
.
assembleResidual
(
_not_ghost
);
}
AKANTU_DEBUG_INFO
(
"Wait distant stresses"
);
// finalize communications
synch_registry
->
waitEndSynchronize
(
_gst_smm_stress
);
AKANTU_DEBUG_INFO
(
"Assemble residual for ghost elements"
);
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
mat
.
assembleResidual
(
_ghost
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
computeStresses
()
{
if
(
isExplicit
())
{
// start synchronization
synch_registry
->
asynchronousSynchronize
(
_gst_smm_uv
);
synch_registry
->
waitEndSynchronize
(
_gst_smm_uv
);
// compute stresses on all local elements for each materials
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
mat
.
computeAllStresses
(
_not_ghost
);
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* Computation of the non local part */
synch_registry
->
asynchronousSynchronize
(
_gst_mnl_for_average
);
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
mat
.
computeAllNonLocalStresses
(
_not_ghost
);
}
synch_registry
->
waitEndSynchronize
(
_gst_mnl_for_average
);
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
mat
.
computeAllNonLocalStresses
(
_ghost
);
}
#endif
}
else
{
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
mat
.
computeAllStressesFromTangentModuli
(
_not_ghost
);
}
}
}
/* -------------------------------------------------------------------------- */
/* Implicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
// /**
// * Initialize the solver and create the sparse matrices needed.
// *
// */
// void SolidMechanicsModel::initSolver(__attribute__((unused))
// SolverOptions & options) {
// UInt nb_global_nodes = mesh.getNbGlobalNodes();
// jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian", _symmetric));
// // jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension);
// if (!isExplicit()) {
// delete stiffness_matrix;
// std::stringstream sstr_sti;
// sstr_sti << id << ":stiffness_matrix";
// stiffness_matrix = &(this->getDOFManager().getNewMatrix("stiffness", _symmetric));
// }
// if (solver) solver->initialize(options);
// }
// /* -------------------------------------------------------------------------- */
// void SolidMechanicsModel::initJacobianMatrix() {
// // @todo make it more flexible: this is an ugly patch to treat the case of non
// // fix profile of the K matrix
// delete jacobian_matrix;
// jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian", "stiffness"));
// std::stringstream sstr_solv; sstr_solv << id << ":solver";
// delete solver;
// solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());
// if(solver)
// solver->initialize(_solver_no_options);
// }
/* -------------------------------------------------------------------------- */
/**
* Initialize the implicit solver, either for dynamic or static cases,
*
* @param dynamic
*/
void
SolidMechanicsModel
::
initImplicit
(
bool
dynamic
)
{
AKANTU_DEBUG_IN
();
method
=
dynamic
?
_implicit_dynamic
:
_static
;
if
(
!
increment
)
setIncrementFlagOn
();
initSolver
();
if
(
method
==
_implicit_dynamic
)
{
if
(
integrator
)
delete
integrator
;
integrator
=
new
TrapezoidalRule2
();
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
assembleStiffnessMatrix
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_INFO
(
"Assemble the new stiffness matrix."
);
// call compute stiffness matrix on each local elements
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
(
*
mat_it
)
->
assembleStiffnessMatrix
(
_not_ghost
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
SparseMatrix
&
SolidMechanicsModel
::
initVelocityDampingMatrix
()
{
if
(
!
velocity_damping_matrix
)
velocity_damping_matrix
=
&
(
this
->
getDOFManager
().
getNewMatrix
(
"velocity_damping"
,
"jacobian"
));
return
*
velocity_damping_matrix
;
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
implicitPred
()
{
AKANTU_DEBUG_IN
();
if
(
previous_displacement
)
previous_displacement
->
copy
(
*
displacement
);
if
(
method
==
_implicit_dynamic
)
integrator
->
integrationSchemePred
(
time_step
,
*
displacement
,
*
velocity
,
*
acceleration
,
*
blocked_dofs
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
implicitCorr
()
{
AKANTU_DEBUG_IN
();
if
(
method
==
_implicit_dynamic
)
{
integrator
->
integrationSchemeCorrDispl
(
time_step
,
*
displacement
,
*
velocity
,
*
acceleration
,
*
blocked_dofs
,
*
increment
);
}
else
{
UInt
nb_nodes
=
displacement
->
getSize
();
UInt
nb_degree_of_freedom
=
displacement
->
getNbComponent
()
*
nb_nodes
;
Real
*
incr_val
=
increment
->
storage
();
Real
*
disp_val
=
displacement
->
storage
();
bool
*
boun_val
=
blocked_dofs
->
storage
();
for
(
UInt
j
=
0
;
j
<
nb_degree_of_freedom
;
++
j
,
++
disp_val
,
++
incr_val
,
++
boun_val
){
*
incr_val
*=
(
1.
-
*
boun_val
);
*
disp_val
+=
*
incr_val
;
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
updateIncrement
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_ASSERT
(
previous_displacement
,
"The previous displacement has to be initialized."
<<
" Are you working with Finite or Ineslactic deformations?"
);
UInt
nb_nodes
=
displacement
->
getSize
();
UInt
nb_degree_of_freedom
=
displacement
->
getNbComponent
()
*
nb_nodes
;
Real
*
incr_val
=
increment
->
storage
();
Real
*
disp_val
=
displacement
->
storage
();
Real
*
prev_disp_val
=
previous_displacement
->
storage
();
for
(
UInt
j
=
0
;
j
<
nb_degree_of_freedom
;
++
j
,
++
disp_val
,
++
incr_val
,
++
prev_disp_val
)
*
incr_val
=
(
*
disp_val
-
*
prev_disp_val
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
updatePreviousDisplacement
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_ASSERT
(
previous_displacement
,
"The previous displacement has to be initialized."
<<
" Are you working with Finite or Ineslactic deformations?"
);
previous_displacement
->
copy
(
*
displacement
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
/* Information */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
synchronizeBoundaries
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_ASSERT
(
synch_registry
,
"Synchronizer registry was not initialized."
<<
" Did you call initParallel?"
);
synch_registry
->
synchronize
(
_gst_smm_boundary
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
synchronizeResidual
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_ASSERT
(
synch_registry
,
"Synchronizer registry was not initialized."
<<
" Did you call initPBC?"
);
synch_registry
->
synchronize
(
_gst_smm_res
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
setIncrementFlagOn
()
{
AKANTU_DEBUG_IN
();
if
(
!
increment
)
{
UInt
nb_nodes
=
mesh
.
getNbNodes
();
std
::
stringstream
sstr_inc
;
sstr_inc
<<
id
<<
":increment"
;
increment
=
&
(
alloc
<
Real
>
(
sstr_inc
.
str
(),
nb_nodes
,
spatial_dimension
,
0.
));
}
increment_flag
=
true
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
Real
SolidMechanicsModel
::
getStableTimeStep
()
{
AKANTU_DEBUG_IN
();
Real
min_dt
=
getStableTimeStep
(
_not_ghost
);
/// reduction min over all processors
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
min_dt
,
1
,
_so_min
);
AKANTU_DEBUG_OUT
();
return
min_dt
;
}
/* -------------------------------------------------------------------------- */
Real
SolidMechanicsModel
::
getStableTimeStep
(
const
GhostType
&
ghost_type
)
{
AKANTU_DEBUG_IN
();
Material
**
mat_val
=
&
(
materials
.
at
(
0
));
Real
min_dt
=
std
::
numeric_limits
<
Real
>::
max
();
updateCurrentPosition
();
Element
elem
;
elem
.
ghost_type
=
ghost_type
;
elem
.
kind
=
_ek_regular
;
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
,
ghost_type
);
Mesh
::
type_iterator
end
=
mesh
.
lastType
(
spatial_dimension
,
ghost_type
);
for
(;
it
!=
end
;
++
it
)
{
elem
.
type
=
*
it
;
UInt
nb_nodes_per_element
=
mesh
.
getNbNodesPerElement
(
*
it
);
UInt
nb_element
=
mesh
.
getNbElement
(
*
it
);
Array
<
UInt
>::
const_scalar_iterator
mat_indexes
=
material_index
(
*
it
,
ghost_type
).
begin
();
Array
<
UInt
>::
const_scalar_iterator
mat_loc_num
=
material_local_numbering
(
*
it
,
ghost_type
).
begin
();
Array
<
Real
>
X
(
0
,
nb_nodes_per_element
*
spatial_dimension
);
FEEngine
::
extractNodalToElementField
(
mesh
,
*
current_position
,
X
,
*
it
,
_not_ghost
);
Array
<
Real
>::
matrix_iterator
X_el
=
X
.
begin
(
spatial_dimension
,
nb_nodes_per_element
);
for
(
UInt
el
=
0
;
el
<
nb_element
;
++
el
,
++
X_el
,
++
mat_indexes
,
++
mat_loc_num
)
{
elem
.
element
=
*
mat_loc_num
;
Real
el_h
=
getFEEngine
().
getElementInradius
(
*
X_el
,
*
it
);
Real
el_c
=
mat_val
[
*
mat_indexes
]
->
getCelerity
(
elem
);
Real
el_dt
=
el_h
/
el_c
;
min_dt
=
std
::
min
(
min_dt
,
el_dt
);
}
}
AKANTU_DEBUG_OUT
();
return
min_dt
;
}
/* -------------------------------------------------------------------------- */
Real
SolidMechanicsModel
::
getKineticEnergy
()
{
AKANTU_DEBUG_IN
();
if
(
!
mass
&&
!
mass_matrix
)
AKANTU_DEBUG_ERROR
(
"No function called to assemble the mass matrix."
);
Real
ekin
=
0.
;
UInt
nb_nodes
=
mesh
.
getNbNodes
();
Real
*
vel_val
=
velocity
->
storage
();
Real
*
mass_val
=
mass
->
storage
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
Real
mv2
=
0
;
bool
is_local_node
=
mesh
.
isLocalOrMasterNode
(
n
);
bool
is_not_pbc_slave_node
=
!
isPBCSlaveNode
(
n
);
bool
count_node
=
is_local_node
&&
is_not_pbc_slave_node
;
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
++
i
)
{
if
(
count_node
)
mv2
+=
*
vel_val
*
*
vel_val
*
*
mass_val
;
vel_val
++
;
mass_val
++
;
}
ekin
+=
mv2
;
}
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
ekin
,
1
,
_so_sum
);
AKANTU_DEBUG_OUT
();
return
ekin
*
.5
;
}
/* -------------------------------------------------------------------------- */
Real
SolidMechanicsModel
::
getKineticEnergy
(
const
ElementType
&
type
,
UInt
index
)
{
AKANTU_DEBUG_IN
();
UInt
nb_quadrature_points
=
getFEEngine
().
getNbQuadraturePoints
(
type
);
Array
<
Real
>
vel_on_quad
(
nb_quadrature_points
,
spatial_dimension
);
Array
<
UInt
>
filter_element
(
1
,
1
,
index
);
getFEEngine
().
interpolateOnQuadraturePoints
(
*
velocity
,
vel_on_quad
,
spatial_dimension
,
type
,
_not_ghost
,
filter_element
);
Array
<
Real
>::
vector_iterator
vit
=
vel_on_quad
.
begin
(
spatial_dimension
);
Array
<
Real
>::
vector_iterator
vend
=
vel_on_quad
.
end
(
spatial_dimension
);
Vector
<
Real
>
rho_v2
(
nb_quadrature_points
);
Real
rho
=
materials
[
material_index
(
type
)(
index
)]
->
getRho
();
for
(
UInt
q
=
0
;
vit
!=
vend
;
++
vit
,
++
q
)
{
rho_v2
(
q
)
=
rho
*
vit
->
dot
(
*
vit
);
}
AKANTU_DEBUG_OUT
();
return
.5
*
getFEEngine
().
integrate
(
rho_v2
,
type
,
index
);
}
/* -------------------------------------------------------------------------- */
Real
SolidMechanicsModel
::
getExternalWork
()
{
AKANTU_DEBUG_IN
();
Real
*
velo
=
velocity
->
storage
();
Real
*
forc
=
force
->
storage
();
Real
*
resi
=
residual
->
storage
();
bool
*
boun
=
blocked_dofs
->
storage
();
Real
work
=
0.
;
UInt
nb_nodes
=
mesh
.
getNbNodes
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
bool
is_local_node
=
mesh
.
isLocalOrMasterNode
(
n
);
bool
is_not_pbc_slave_node
=
!
isPBCSlaveNode
(
n
);
bool
count_node
=
is_local_node
&&
is_not_pbc_slave_node
;
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
++
i
)
{
if
(
count_node
)
{
if
(
*
boun
)
work
-=
*
resi
*
*
velo
*
time_step
;
else
work
+=
*
forc
*
*
velo
*
time_step
;
}
++
velo
;
++
forc
;
++
resi
;
++
boun
;
}
}
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
work
,
1
,
_so_sum
);
AKANTU_DEBUG_OUT
();
return
work
;
}
/* -------------------------------------------------------------------------- */
Real
SolidMechanicsModel
::
getEnergy
(
const
std
::
string
&
energy_id
)
{
AKANTU_DEBUG_IN
();
if
(
energy_id
==
"kinetic"
)
{
return
getKineticEnergy
();
}
else
if
(
energy_id
==
"external work"
){
return
getExternalWork
();
}
Real
energy
=
0.
;
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
energy
+=
(
*
mat_it
)
->
getEnergy
(
energy_id
);
}
/// reduction sum over all processors
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
energy
,
1
,
_so_sum
);
AKANTU_DEBUG_OUT
();
return
energy
;
}
/* -------------------------------------------------------------------------- */
Real
SolidMechanicsModel
::
getEnergy
(
const
std
::
string
&
energy_id
,
const
ElementType
&
type
,
UInt
index
){
AKANTU_DEBUG_IN
();
if
(
energy_id
==
"kinetic"
)
{
return
getKineticEnergy
(
type
,
index
);
}
std
::
vector
<
Material
*>::
iterator
mat_it
;
UInt
mat_index
=
this
->
material_index
(
type
,
_not_ghost
)(
index
);
UInt
mat_loc_num
=
this
->
material_local_numbering
(
type
,
_not_ghost
)(
index
);
Real
energy
=
this
->
materials
[
mat_index
]
->
getEnergy
(
energy_id
,
type
,
mat_loc_num
);
AKANTU_DEBUG_OUT
();
return
energy
;
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
onElementsAdded
(
const
Array
<
Element
>
&
element_list
,
const
NewElementsEvent
&
event
)
{
AKANTU_DEBUG_IN
();
this
->
getFEEngine
().
initShapeFunctions
(
_not_ghost
);
this
->
getFEEngine
().
initShapeFunctions
(
_ghost
);
for
(
UInt
g
=
_not_ghost
;
g
<=
_ghost
;
++
g
)
{
GhostType
gt
=
(
GhostType
)
g
;
Mesh
::
type_iterator
it
=
this
->
mesh
.
firstType
(
spatial_dimension
,
gt
,
_ek_not_defined
);
Mesh
::
type_iterator
end
=
this
->
mesh
.
lastType
(
spatial_dimension
,
gt
,
_ek_not_defined
);
for
(;
it
!=
end
;
++
it
)
{
UInt
nb_element
=
this
->
mesh
.
getNbElement
(
*
it
,
gt
);
if
(
!
material_index
.
exists
(
*
it
,
gt
))
{
this
->
material_index
.
alloc
(
nb_element
,
1
,
*
it
,
gt
);
this
->
material_local_numbering
.
alloc
(
nb_element
,
1
,
*
it
,
gt
);
}
else
{
this
->
material_index
(
*
it
,
gt
).
resize
(
nb_element
);
this
->
material_local_numbering
(
*
it
,
gt
).
resize
(
nb_element
);
}
}
}
Array
<
Element
>::
const_iterator
<
Element
>
it
=
element_list
.
begin
();
Array
<
Element
>::
const_iterator
<
Element
>
end
=
element_list
.
end
();
ElementTypeMapArray
<
UInt
>
filter
(
"new_element_filter"
,
this
->
getID
());
for
(
UInt
el
=
0
;
it
!=
end
;
++
it
,
++
el
)
{
const
Element
&
elem
=
*
it
;
if
(
!
filter
.
exists
(
elem
.
type
,
elem
.
ghost_type
))
filter
.
alloc
(
0
,
1
,
elem
.
type
,
elem
.
ghost_type
);
filter
(
elem
.
type
,
elem
.
ghost_type
).
push_back
(
elem
.
element
);
}
this
->
assignMaterialToElements
(
&
filter
);
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
(
*
mat_it
)
->
onElementsAdded
(
element_list
,
event
);
}
if
(
method
==
_explicit_lumped_mass
)
this
->
assembleMassLumped
();
if
(
method
!=
_explicit_lumped_mass
)
{
this
->
initSolver
();
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
onElementsRemoved
(
__attribute__
((
unused
))
const
Array
<
Element
>
&
element_list
,
const
ElementTypeMapArray
<
UInt
>
&
new_numbering
,
const
RemovedElementsEvent
&
event
)
{
this
->
getFEEngine
().
initShapeFunctions
(
_not_ghost
);
this
->
getFEEngine
().
initShapeFunctions
(
_ghost
);
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
(
*
mat_it
)
->
onElementsRemoved
(
element_list
,
new_numbering
,
event
);
}
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
onNodesAdded
(
const
Array
<
UInt
>
&
nodes_list
,
__attribute__
((
unused
))
const
NewNodesEvent
&
event
)
{
AKANTU_DEBUG_IN
();
UInt
nb_nodes
=
mesh
.
getNbNodes
();
if
(
displacement
)
displacement
->
resize
(
nb_nodes
);
if
(
mass
)
mass
->
resize
(
nb_nodes
);
if
(
velocity
)
velocity
->
resize
(
nb_nodes
);
if
(
acceleration
)
acceleration
->
resize
(
nb_nodes
);
if
(
force
)
force
->
resize
(
nb_nodes
);
if
(
residual
)
residual
->
resize
(
nb_nodes
);
if
(
blocked_dofs
)
blocked_dofs
->
resize
(
nb_nodes
);
if
(
previous_displacement
)
previous_displacement
->
resize
(
nb_nodes
);
if
(
increment_acceleration
)
increment_acceleration
->
resize
(
nb_nodes
);
if
(
increment
)
increment
->
resize
(
nb_nodes
);
if
(
current_position
)
current_position
->
resize
(
nb_nodes
);
delete
dof_synchronizer
;
dof_synchronizer
=
new
DOFSynchronizer
(
mesh
,
spatial_dimension
);
dof_synchronizer
->
initLocalDOFEquationNumbers
();
dof_synchronizer
->
initGlobalDOFEquationNumbers
();
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
(
*
mat_it
)
->
onNodesAdded
(
nodes_list
,
event
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
onNodesRemoved
(
__attribute__
((
unused
))
const
Array
<
UInt
>
&
element_list
,
const
Array
<
UInt
>
&
new_numbering
,
__attribute__
((
unused
))
const
RemovedNodesEvent
&
event
)
{
if
(
displacement
)
mesh
.
removeNodesFromArray
(
*
displacement
,
new_numbering
);
if
(
mass
)
mesh
.
removeNodesFromArray
(
*
mass
,
new_numbering
);
if
(
velocity
)
mesh
.
removeNodesFromArray
(
*
velocity
,
new_numbering
);
if
(
acceleration
)
mesh
.
removeNodesFromArray
(
*
acceleration
,
new_numbering
);
if
(
force
)
mesh
.
removeNodesFromArray
(
*
force
,
new_numbering
);
if
(
residual
)
mesh
.
removeNodesFromArray
(
*
residual
,
new_numbering
);
if
(
blocked_dofs
)
mesh
.
removeNodesFromArray
(
*
blocked_dofs
,
new_numbering
);
if
(
increment_acceleration
)
mesh
.
removeNodesFromArray
(
*
increment_acceleration
,
new_numbering
);
if
(
increment
)
mesh
.
removeNodesFromArray
(
*
increment
,
new_numbering
);
delete
dof_synchronizer
;
dof_synchronizer
=
new
DOFSynchronizer
(
mesh
,
spatial_dimension
);
dof_synchronizer
->
initLocalDOFEquationNumbers
();
dof_synchronizer
->
initGlobalDOFEquationNumbers
();
if
(
method
!=
_explicit_lumped_mass
)
{
this
->
initSolver
();
}
}
/* -------------------------------------------------------------------------- */
bool
SolidMechanicsModel
::
isInternal
(
const
std
::
string
&
field_name
,
const
ElementKind
&
element_kind
){
bool
is_internal
=
false
;
/// check if at least one material contains field_id as an internal
for
(
UInt
m
=
0
;
m
<
materials
.
size
()
&&
!
is_internal
;
++
m
)
{
is_internal
|=
materials
[
m
]
->
isInternal
(
field_name
,
element_kind
);
}
return
is_internal
;
}
/* -------------------------------------------------------------------------- */
ElementTypeMap
<
UInt
>
SolidMechanicsModel
::
getInternalDataPerElem
(
const
std
::
string
&
field_name
,
const
ElementKind
&
element_kind
){
if
(
!
(
this
->
isInternal
(
field_name
,
element_kind
)))
AKANTU_EXCEPTION
(
"unknown internal "
<<
field_name
);
for
(
UInt
m
=
0
;
m
<
materials
.
size
()
;
++
m
)
{
if
(
materials
[
m
]
->
isInternal
(
field_name
,
element_kind
))
return
materials
[
m
]
->
getInternalDataPerElem
(
field_name
,
element_kind
);
}
return
ElementTypeMap
<
UInt
>
();
}
/* -------------------------------------------------------------------------- */
ElementTypeMapArray
<
Real
>
&
SolidMechanicsModel
::
flattenInternal
(
const
std
::
string
&
field_name
,
const
ElementKind
&
kind
,
const
GhostType
ghost_type
){
std
::
pair
<
std
::
string
,
ElementKind
>
key
(
field_name
,
kind
);
if
(
this
->
registered_internals
.
count
(
key
)
==
0
){
this
->
registered_internals
[
key
]
=
new
ElementTypeMapArray
<
Real
>
(
field_name
,
this
->
id
);
}
ElementTypeMapArray
<
Real
>
*
internal_flat
=
this
->
registered_internals
[
key
];
typedef
ElementTypeMapArray
<
Real
>::
type_iterator
iterator
;
iterator
tit
=
internal_flat
->
firstType
(
spatial_dimension
,
ghost_type
,
kind
);
iterator
end
=
internal_flat
->
lastType
(
spatial_dimension
,
ghost_type
,
kind
);
for
(;
tit
!=
end
;
++
tit
)
{
ElementType
type
=
*
tit
;
(
*
internal_flat
)(
type
,
ghost_type
).
clear
();
}
for
(
UInt
m
=
0
;
m
<
materials
.
size
();
++
m
)
{
if
(
materials
[
m
]
->
isInternal
(
field_name
,
kind
))
materials
[
m
]
->
flattenInternal
(
field_name
,
*
internal_flat
,
ghost_type
,
kind
);
}
return
*
internal_flat
;
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
flattenAllRegisteredInternals
(
const
ElementKind
&
kind
){
std
::
map
<
std
::
pair
<
std
::
string
,
ElementKind
>
,
ElementTypeMapArray
<
Real
>
*>::
iterator
it
=
this
->
registered_internals
.
begin
();
std
::
map
<
std
::
pair
<
std
::
string
,
ElementKind
>
,
ElementTypeMapArray
<
Real
>
*>::
iterator
end
=
this
->
registered_internals
.
end
();
while
(
it
!=
end
){
if
(
kind
==
it
->
first
.
second
)
this
->
flattenInternal
(
it
->
first
.
first
,
kind
);
++
it
;
}
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
onDump
(){
this
->
flattenAllRegisteredInternals
(
_ek_regular
);
}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
dumper
::
Field
*
SolidMechanicsModel
::
createElementalField
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
,
const
UInt
&
spatial_dimension
,
const
ElementKind
&
kind
)
{
dumper
::
Field
*
field
=
NULL
;
if
(
field_name
==
"partitions"
)
field
=
mesh
.
createElementalField
<
UInt
,
dumper
::
ElementPartitionField
>
(
mesh
.
getConnectivities
(),
group_name
,
spatial_dimension
,
kind
);
else
if
(
field_name
==
"material_index"
)
field
=
mesh
.
createElementalField
<
UInt
,
Vector
,
dumper
::
ElementalField
>
(
material_index
,
group_name
,
spatial_dimension
,
kind
);
else
{
// this copy of field_name is used to compute derivated data such as
// strain and von mises stress that are based on grad_u and stress
std
::
string
field_name_copy
(
field_name
);
if
(
field_name
==
"strain"
||
field_name
==
"Green strain"
||
field_name
==
"principal strain"
||
field_name
==
"principal Green strain"
)
field_name_copy
=
"grad_u"
;
else
if
(
field_name
==
"Von Mises stress"
)
field_name_copy
=
"stress"
;
bool
is_internal
=
this
->
isInternal
(
field_name_copy
,
kind
);
if
(
is_internal
)
{
ElementTypeMap
<
UInt
>
nb_data_per_elem
=
this
->
getInternalDataPerElem
(
field_name_copy
,
kind
);
ElementTypeMapArray
<
Real
>
&
internal_flat
=
this
->
flattenInternal
(
field_name_copy
,
kind
);
field
=
mesh
.
createElementalField
<
Real
,
dumper
::
InternalMaterialField
>
(
internal_flat
,
group_name
,
spatial_dimension
,
kind
,
nb_data_per_elem
);
if
(
field_name
==
"strain"
){
dumper
::
ComputeStrain
<
false
>
*
foo
=
new
dumper
::
ComputeStrain
<
false
>
(
*
this
);
field
=
dumper
::
FieldComputeProxy
::
createFieldCompute
(
field
,
*
foo
);
}
else
if
(
field_name
==
"Von Mises stress"
)
{
dumper
::
ComputeVonMisesStress
*
foo
=
new
dumper
::
ComputeVonMisesStress
(
*
this
);
field
=
dumper
::
FieldComputeProxy
::
createFieldCompute
(
field
,
*
foo
);
}
else
if
(
field_name
==
"Green strain"
)
{
dumper
::
ComputeStrain
<
true
>
*
foo
=
new
dumper
::
ComputeStrain
<
true
>
(
*
this
);
field
=
dumper
::
FieldComputeProxy
::
createFieldCompute
(
field
,
*
foo
);
}
else
if
(
field_name
==
"principal strain"
)
{
dumper
::
ComputePrincipalStrain
<
false
>
*
foo
=
new
dumper
::
ComputePrincipalStrain
<
false
>
(
*
this
);
field
=
dumper
::
FieldComputeProxy
::
createFieldCompute
(
field
,
*
foo
);
}
else
if
(
field_name
==
"principal Green strain"
)
{
dumper
::
ComputePrincipalStrain
<
true
>
*
foo
=
new
dumper
::
ComputePrincipalStrain
<
true
>
(
*
this
);
field
=
dumper
::
FieldComputeProxy
::
createFieldCompute
(
field
,
*
foo
);
}
//treat the paddings
if
(
padding_flag
){
if
(
field_name
==
"stress"
){
if
(
spatial_dimension
==
2
)
{
dumper
::
StressPadder
<
2
>
*
foo
=
new
dumper
::
StressPadder
<
2
>
(
*
this
);
field
=
dumper
::
FieldComputeProxy
::
createFieldCompute
(
field
,
*
foo
);
}
}
else
if
(
field_name
==
"strain"
||
field_name
==
"Green strain"
){
if
(
spatial_dimension
==
2
)
{
dumper
::
StrainPadder
<
2
>
*
foo
=
new
dumper
::
StrainPadder
<
2
>
(
*
this
);
field
=
dumper
::
FieldComputeProxy
::
createFieldCompute
(
field
,
*
foo
);
}
}
}
// homogenize the field
dumper
::
ComputeFunctorInterface
*
foo
=
dumper
::
HomogenizerProxy
::
createHomogenizer
(
*
field
);
field
=
dumper
::
FieldComputeProxy
::
createFieldCompute
(
field
,
*
foo
);
}
}
return
field
;
}
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
SolidMechanicsModel
::
createNodalFieldReal
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
)
{
std
::
map
<
std
::
string
,
Array
<
Real
>*
>
real_nodal_fields
;
real_nodal_fields
[
"displacement"
]
=
displacement
;
real_nodal_fields
[
"mass"
]
=
mass
;
real_nodal_fields
[
"velocity"
]
=
velocity
;
real_nodal_fields
[
"acceleration"
]
=
acceleration
;
real_nodal_fields
[
"force"
]
=
force
;
real_nodal_fields
[
"residual"
]
=
residual
;
real_nodal_fields
[
"increment"
]
=
increment
;
dumper
::
Field
*
field
=
NULL
;
if
(
padding_flag
)
field
=
mesh
.
createNodalField
(
real_nodal_fields
[
field_name
],
group_name
,
3
);
else
field
=
mesh
.
createNodalField
(
real_nodal_fields
[
field_name
],
group_name
);
return
field
;
}
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
SolidMechanicsModel
::
createNodalFieldBool
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
)
{
std
::
map
<
std
::
string
,
Array
<
bool
>*
>
uint_nodal_fields
;
uint_nodal_fields
[
"blocked_dofs"
]
=
blocked_dofs
;
dumper
::
Field
*
field
=
NULL
;
field
=
mesh
.
createNodalField
(
uint_nodal_fields
[
field_name
],
group_name
);
return
field
;
}
/* -------------------------------------------------------------------------- */
#else
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
SolidMechanicsModel
::
createElementalField
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
,
const
UInt
&
spatial_dimension
,
const
ElementKind
&
kind
){
return
NULL
;
}
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
SolidMechanicsModel
::
createNodalFieldReal
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
)
{
return
NULL
;
}
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
SolidMechanicsModel
::
createNodalFieldBool
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
)
{
return
NULL
;
}
#endif
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
dump
(
const
std
::
string
&
dumper_name
)
{
this
->
onDump
();
EventManager
::
sendEvent
(
SolidMechanicsModelEvent
::
BeforeDumpEvent
());
mesh
.
dump
(
dumper_name
);
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
dump
(
const
std
::
string
&
dumper_name
,
UInt
step
)
{
this
->
onDump
();
EventManager
::
sendEvent
(
SolidMechanicsModelEvent
::
BeforeDumpEvent
());
mesh
.
dump
(
dumper_name
,
step
);
}
/* ------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
dump
(
const
std
::
string
&
dumper_name
,
Real
time
,
UInt
step
)
{
this
->
onDump
();
EventManager
::
sendEvent
(
SolidMechanicsModelEvent
::
BeforeDumpEvent
());
mesh
.
dump
(
dumper_name
,
time
,
step
);
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
dump
()
{
this
->
onDump
();
EventManager
::
sendEvent
(
SolidMechanicsModelEvent
::
BeforeDumpEvent
());
mesh
.
dump
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
dump
(
UInt
step
)
{
this
->
onDump
();
EventManager
::
sendEvent
(
SolidMechanicsModelEvent
::
BeforeDumpEvent
());
mesh
.
dump
(
step
);
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
dump
(
Real
time
,
UInt
step
)
{
this
->
onDump
();
EventManager
::
sendEvent
(
SolidMechanicsModelEvent
::
BeforeDumpEvent
());
mesh
.
dump
(
time
,
step
);
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
computeCauchyStresses
()
{
AKANTU_DEBUG_IN
();
// call compute stiffness matrix on each local elements
std
::
vector
<
Material
*>::
iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
Material
&
mat
=
**
mat_it
;
if
(
mat
.
isFiniteDeformation
())
mat
.
computeAllCauchyStresses
(
_not_ghost
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
saveStressAndStrainBeforeDamage
()
{
EventManager
::
sendEvent
(
SolidMechanicsModelEvent
::
BeginningOfDamageIterationEvent
());
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
updateEnergiesAfterDamage
()
{
EventManager
::
sendEvent
(
SolidMechanicsModelEvent
::
AfterDamageEvent
());
}
/* -------------------------------------------------------------------------- */
void
SolidMechanicsModel
::
printself
(
std
::
ostream
&
stream
,
int
indent
)
const
{
std
::
string
space
;
for
(
Int
i
=
0
;
i
<
indent
;
i
++
,
space
+=
AKANTU_INDENT
);
stream
<<
space
<<
"Solid Mechanics Model ["
<<
std
::
endl
;
stream
<<
space
<<
" + id : "
<<
id
<<
std
::
endl
;
stream
<<
space
<<
" + spatial dimension : "
<<
spatial_dimension
<<
std
::
endl
;
stream
<<
space
<<
" + fem ["
<<
std
::
endl
;
getFEEngine
().
printself
(
stream
,
indent
+
2
);
stream
<<
space
<<
AKANTU_INDENT
<<
"]"
<<
std
::
endl
;
stream
<<
space
<<
" + nodals information ["
<<
std
::
endl
;
displacement
->
printself
(
stream
,
indent
+
2
);
mass
->
printself
(
stream
,
indent
+
2
);
velocity
->
printself
(
stream
,
indent
+
2
);
acceleration
->
printself
(
stream
,
indent
+
2
);
force
->
printself
(
stream
,
indent
+
2
);
residual
->
printself
(
stream
,
indent
+
2
);
blocked_dofs
->
printself
(
stream
,
indent
+
2
);
stream
<<
space
<<
AKANTU_INDENT
<<
"]"
<<
std
::
endl
;
stream
<<
space
<<
" + material information ["
<<
std
::
endl
;
material_index
.
printself
(
stream
,
indent
+
2
);
stream
<<
space
<<
AKANTU_INDENT
<<
"]"
<<
std
::
endl
;
stream
<<
space
<<
" + materials ["
<<
std
::
endl
;
std
::
vector
<
Material
*>::
const_iterator
mat_it
;
for
(
mat_it
=
materials
.
begin
();
mat_it
!=
materials
.
end
();
++
mat_it
)
{
const
Material
&
mat
=
*
(
*
mat_it
);
mat
.
printself
(
stream
,
indent
+
1
);
}
stream
<<
space
<<
AKANTU_INDENT
<<
"]"
<<
std
::
endl
;
stream
<<
space
<<
"]"
<<
std
::
endl
;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
Event Timeline
Log In to Comment