Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82324450
domain_akantu_dynamic.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
Tue, Sep 10, 18:47
Size
10 KB
Mime Type
text/x-c++
Expires
Thu, Sep 12, 18:47 (2 d)
Engine
blob
Format
Raw Data
Handle
20672210
Attached To
rLIBMULTISCALE LibMultiScale
domain_akantu_dynamic.cc
View Options
/**
* @file domain_akantu_dynamic.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Tue Jul 22 14:47:56 2014
*
* @brief This is the model wrapping Akantu in its dynamic solve version
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale 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.
*
* LibMultiScale 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 LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#define TIMER
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "domain_akantu_dynamic.hh"
#include "lib_geometry.hh"
#include "math_tools.hh"
#include "units.hh"
/* -------------------------------------------------------------------------- */
#include <mesh_io.hh>
#include <solid_mechanics_model.hh>
#ifdef AKANTU_HEAT_TRANSFER
#include <heat_transfer_model.hh>
#endif
#include <material.hh>
#include <material_elastic.hh>
/* -------------------------------------------------------------------------- */
#include <iomanip>
#include <sstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
>
DomainAkantuDynamic
<
Dim
>::~
DomainAkantuDynamic
(){
#ifdef AKANTU_HEAT_TRANSFER
delete
(
this
->
heat_transfer_model
);
this
->
heat_transfer_model
=
NULL
;
delete
(
this
->
temperature
);
this
->
temperature
=
NULL
;
delete
(
this
->
temperature_rate
);
this
->
temperature_rate
=
NULL
;
delete
(
this
->
heat_flux
);
this
->
heat_flux
=
NULL
;
#endif
}
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
>
DomainAkantuDynamic
<
Dim
>::
DomainAkantuDynamic
(
DomainID
ID
,
CommGroup
GID
)
:
DomainAkantu
<
Dim
>
(
ID
,
GID
),
heat_transfer_model
(
NULL
),
heat_transfer_flag
(
false
),
temperature
(
NULL
),
temperature_rate
(
NULL
),
heat_flux
(
NULL
),
heat_boundary_geom
(
invalidGeom
)
{}
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
>
void
DomainAkantuDynamic
<
Dim
>::
init
(){
if
(
this
->
material_filename
!=
""
)
this
->
akantu_parser
.
parse
(
this
->
material_filename
);
this
->
model
->
setParser
(
this
->
akantu_parser
);
this
->
model
->
initFull
(
akantu
::
SolidMechanicsModelOptions
(
akantu
::
_explicit_lumped_mass
,
true
));
this
->
initMaterial
();
akantu
::
Array
<
bool
>
&
bound
=
this
->
model
->
getBlockedDOFs
();
UInt
full_counter
=
0
;
for
(
UInt
i
=
0
;
i
<
bound
.
getSize
()
;
++
i
)
{
if
(
bound
(
i
))
{
full_counter
++
;
}
}
this
->
elems
.
setSolidMechanicsModel
(
*
this
->
model
);
this
->
nodes
.
setSolidMechanicsModel
(
*
this
->
model
);
Real
time_step_suggested
=
this
->
model
->
getStableTimeStep
()
/
sqrt
(
code_unit_system
.
e_m2dd_tt
);
DUMP
(
"suggested timestep is for solid mechanics model "
<<
time_step_suggested
,
DBG_MESSAGE
);
if
(
this
->
timeStep
>
time_step_suggested
)
LM_FATAL
(
"provided timestep is larger than the critical timestep: please adapt it"
);
this
->
model
->
setTimeStep
(
this
->
timeStep
);
}
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
>
void
DomainAkantuDynamic
<
Dim
>::
initHeatTransfer
(){
#ifndef AKANTU_HEAT_TRANSFER
LM_FATAL
(
"Akantu was not compiled using the heat transfer model"
);
#else
heat_transfer_model
=
new
akantu
::
HeatTransferModel
(
*
this
->
mesh
,
Dim
,
this
->
getID
()
+
"heat"
);
//initialize PBC
bool
pbc_activated
=
false
;
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
pbc_activated
|=
this
->
pbc
[
i
];
}
if
(
pbc_activated
)
{
heat_transfer_model
->
setPBC
(
this
->
pbc
[
0
],
this
->
pbc
[
1
],
this
->
pbc
[
2
]);
std
::
map
<
UInt
,
UInt
>
&
akantu_pairs
=
heat_transfer_model
->
getPBCPairs
();
std
::
map
<
UInt
,
UInt
>::
iterator
it
=
akantu_pairs
.
begin
();
std
::
map
<
UInt
,
UInt
>::
iterator
end
=
akantu_pairs
.
end
();
while
(
it
!=
end
)
{
this
->
pbc_pairs
.
push_back
((
*
it
));
++
it
;
}
}
//initialize everything
if
(
heat_transfer_material_filename
!=
""
)
heat_transfer_model
->
initFull
();
else
LM_FATAL
(
"Heat transfer material file not provided."
);
this
->
elems
.
setHeatTransferModel
(
*
heat_transfer_model
);
this
->
nodes
.
setHeatTransferModel
(
*
heat_transfer_model
);
this
->
elems
.
setHeatTransferFlag
(
heat_transfer_flag
);
this
->
nodes
.
setHeatTransferFlag
(
heat_transfer_flag
);
//assemble the lumped capacity
heat_transfer_model
->
assembleCapacityLumped
();
//get stable time step
akantu
::
Real
suggested_time_step_heat_transfer
=
heat_transfer_model
->
getStableTimeStep
()
*
0.8
;
DUMP
(
"suggested timestep for heat transfer model is "
<<
suggested_time_step_heat_transfer
,
DBG_MESSAGE
);
if
(
this
->
timeStep
>
suggested_time_step_heat_transfer
)
LM_FATAL
(
"provided timestep is larger than the critical timestep for the heat transfer model: please adapt it"
);
heat_transfer_model
->
setTimeStep
(
this
->
timeStep
);
/* -------------------- */
// initialize the vectors
/* -------------------- */
temperature
=
new
VecAkantu
(
heat_transfer_model
->
getTemperature
());
heat_flux
=
new
VecAkantu
(
heat_transfer_model
->
getExternalHeatRate
());
temperature_rate
=
new
VecAkantu
(
heat_transfer_model
->
getTemperatureRate
());
heat_transfer_model
->
updateResidual
();
akantu
::
Array
<
bool
>
&
boundary
=
heat_transfer_model
->
getBlockedDOFs
();
UInt
nb_nodes
=
boundary
.
getSize
();
bool
*
bound
=
boundary
.
storage
();
if
(
heat_boundary_geom
!=
invalidGeom
){
Geometry
*
geom_dirichlet
=
GeometryManager
::
getManager
().
getGeometry
(
heat_boundary_geom
);
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
,
++
bound
){
Real
pos
[
Dim
];
this
->
nodes
.
get
(
n
).
getPositions0
(
pos
);
if
(
geom_dirichlet
->
contains
<
Dim
>
(
pos
)){
*
bound
=
true
;
}
}
}
#endif
}
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
>
void
DomainAkantuDynamic
<
Dim
>::
performStep1
(){
STARTTIMER
(
"AkantuStepPosition"
);
DomainAkantu
<
Dim
>::
performStep1
();
UInt
nb_nodes
=
this
->
model
->
getFEEngine
().
getMesh
().
getNbNodes
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
(
*
this
->
velocity
)[
Dim
*
n
+
i
]
+=
code_unit_system
.
ft_m2v
*
.5
*
(
*
this
->
residual
)[
Dim
*
n
+
i
]
/
(
*
this
->
mass
)[
Dim
*
n
+
i
]
*
this
->
timeStep
;
LM_ASSERT
(
!
isnan
((
*
this
->
velocity
)[
Dim
*
n
+
i
])
&&
!
isinf
((
*
this
->
velocity
)[
Dim
*
n
+
i
]),
"problem with velocity on node "
<<
n
<<
" Dim "
<<
i
<<
" "
<<
(
*
this
->
residual
)[
Dim
*
n
+
i
]
<<
" "
<<
(
*
this
->
mass
)[
Dim
*
n
+
i
]);
(
*
this
->
displacement
)[
Dim
*
n
+
i
]
+=
this
->
timeStep
*
(
*
this
->
velocity
)[
Dim
*
n
+
i
];
}
}
STOPTIMER
(
"AkantuStepPosition"
);
#ifdef AKANTU_HEAT_TRANSFER
if
(
heat_transfer_flag
)
heat_transfer_model
->
explicitPred
();
#endif
this
->
mesh_container
.
incRelease
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
>
void
DomainAkantuDynamic
<
Dim
>::
performStep2
(){
STARTTIMER
(
"AkantuStepForce"
);
this
->
model
->
initializeUpdateResidualData
();
if
(
this
->
prestressed_periodicity_flag
)
{
typedef
std
::
vector
<
std
::
pair
<
UInt
,
UInt
>
>::
iterator
pair_it
;
pair_it
pbc_pair
=
this
->
pbc_pairs
.
begin
();
pair_it
end
=
this
->
pbc_pairs
.
end
();
UInt
node_counter
=
0
;
for
(;
pbc_pair
!=
end
;
++
pbc_pair
,
++
node_counter
)
{
for
(
UInt
direction
=
0
;
direction
<
Dim
;
++
direction
)
{
(
*
this
->
displacement
)[
Dim
*
pbc_pair
->
first
+
direction
]
+=
this
->
slave_displacements
[
Dim
*
node_counter
+
direction
];
}
}
}
this
->
model
->
updateResidual
(
false
);
// model->updateAcceleration();
STOPTIMER
(
"AkantuStepForce"
);
#ifdef AKANTU_HEAT_TRANSFER
if
(
heat_transfer_flag
)
heat_transfer_model
->
updateResidual
();
#endif
this
->
mesh_container
.
incRelease
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
Dim
>
void
DomainAkantuDynamic
<
Dim
>::
performStep3
(){
STARTTIMER
(
"AkantuStepVelocity"
);
DomainAkantu
<
Dim
>::
performStep3
();
UInt
nb_nodes
=
this
->
model
->
getFEEngine
().
getMesh
().
getNbNodes
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
for
(
UInt
i
=
0
;
i
<
Dim
;
++
i
)
{
(
*
this
->
velocity
)[
Dim
*
n
+
i
]
+=
code_unit_system
.
ft_m2v
*
.5
*
(
*
this
->
residual
)[
Dim
*
n
+
i
]
/
(
*
this
->
mass
)[
Dim
*
n
+
i
]
*
this
->
timeStep
;
LM_ASSERT
(
!
isnan
((
*
this
->
velocity
)[
Dim
*
n
+
i
])
&&
!
isinf
((
*
this
->
velocity
)[
Dim
*
n
+
i
]),
"problem with velocity on node "
<<
n
<<
" Dim "
<<
i
<<
" "
<<
(
*
this
->
velocity
)[
Dim
*
n
+
i
]
<<
" "
<<
(
*
this
->
residual
)[
Dim
*
n
+
i
]
<<
" "
<<
(
*
this
->
mass
)[
Dim
*
n
+
i
]
<<
" "
<<
this
->
timeStep
);
}
}
STOPTIMER
(
"AkantuStepVelocity"
);
#ifdef AKANTU_HEAT_TRANSFER
if
(
heat_transfer_flag
)
heat_transfer_model
->
explicitCorr
();
#endif
this
->
mesh_container
.
incRelease
();
}
/* -------------------------------------------------------------------------- */
/* LMDESC AKANTU_DYNAMIC
This domain implements the plugin with Akantu Finite Element library.
*/
/* LMHERITANCE domain_akantu */
template
<
UInt
Dim
>
void
DomainAkantuDynamic
<
Dim
>::
declareParams
(){
//DomainAkantu<Dim>::declareParams();
/* LMKEYWORD HEAT_TRANSFER
To switch on the heat transfer model object
*/
this
->
parseTag
(
"HEAT_TRANSFER"
,
this
->
heat_transfer_flag
,
false
);
/* LMKEYWORD HEAT_TRANSFER_MATERIAL_FILENAME
Specify the material file to be loaded for the Akantu HeatTransferModel object
*/
this
->
parseKeyword
(
"HEAT_TRANSFER_MATERIAL_FILENAME"
,
heat_transfer_material_filename
,
""
);
/* LMKEYWORD HEAT_BOUNDARY
Specify the geometry for which a Dirichlet boundary condition should be used for
the HeatTransfer model
*/
this
->
parseKeyword
(
"HEAT_BOUNDARY"
,
heat_boundary_geom
,
invalidGeom
);
}
template
class
DomainAkantuDynamic
<
1
>
;
template
class
DomainAkantuDynamic
<
2
>
;
template
class
DomainAkantuDynamic
<
3
>
;
__END_LIBMULTISCALE__
Event Timeline
Log In to Comment