Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80760995
heat_transfer_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, Sep 2, 11:26
Size
40 KB
Mime Type
text/x-c
Expires
Wed, Sep 4, 11:26 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20402071
Attached To
rAKA akantu
heat_transfer_model.cc
View Options
/**
* @file heat_transfer_model.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Lucas Frerot <lucas.frerot@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Rui Wang <rui.wang@epfl.ch>
*
* @date creation: Sun May 01 2011
* @date last modification: Mon Nov 30 2015
*
* @brief Implementation of HeatTransferModel class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "heat_transfer_model.hh"
#include "group_manager_inline_impl.cc"
#include "dumpable_inline_impl.hh"
#include "aka_math.hh"
#include "aka_common.hh"
#include "fe_engine_template.hh"
#include "mesh.hh"
#include "static_communicator.hh"
#include "parser.hh"
#include "generalized_trapezoidal.hh"
#ifdef AKANTU_USE_MUMPS
#include "solver_mumps.hh"
#endif
#ifdef AKANTU_USE_IOHELPER
# include "dumper_paraview.hh"
# include "dumper_elemental_field.hh"
# include "dumper_element_partition.hh"
# include "dumper_internal_material_field.hh"
#endif
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
const
HeatTransferModelOptions
default_heat_transfer_model_options
(
_explicit_lumped_capacity
);
/* -------------------------------------------------------------------------- */
HeatTransferModel
::
HeatTransferModel
(
Mesh
&
mesh
,
UInt
dim
,
const
ID
&
id
,
const
MemoryID
&
memory_id
)
:
Model
(
mesh
,
dim
,
id
,
memory_id
),
Parsable
(
_st_heat
,
id
),
integrator
(
NULL
),
conductivity_matrix
(
NULL
),
capacity_matrix
(
NULL
),
jacobian_matrix
(
NULL
),
temperature_gradient
(
"temperature_gradient"
,
id
),
temperature_on_qpoints
(
"temperature_on_qpoints"
,
id
),
conductivity_on_qpoints
(
"conductivity_on_qpoints"
,
id
),
k_gradt_on_qpoints
(
"k_gradt_on_qpoints"
,
id
),
int_bt_k_gT
(
"int_bt_k_gT"
,
id
),
bt_k_gT
(
"bt_k_gT"
,
id
),
conductivity
(
spatial_dimension
,
spatial_dimension
),
thermal_energy
(
"thermal_energy"
,
id
),
solver
(
NULL
),
pbc_synch
(
NULL
)
{
AKANTU_DEBUG_IN
();
createSynchronizerRegistry
(
this
);
std
::
stringstream
sstr
;
sstr
<<
id
<<
":fem"
;
registerFEEngineObject
<
MyFEEngineType
>
(
sstr
.
str
(),
mesh
,
spatial_dimension
);
this
->
temperature
=
NULL
;
this
->
residual
=
NULL
;
this
->
blocked_dofs
=
NULL
;
#ifdef AKANTU_USE_IOHELPER
this
->
mesh
.
registerDumper
<
DumperParaview
>
(
"paraview_all"
,
id
,
true
);
this
->
mesh
.
addDumpMesh
(
mesh
,
spatial_dimension
,
_not_ghost
,
_ek_regular
);
#endif
this
->
registerParam
(
"conductivity"
,
conductivity
,
_pat_parsmod
);
this
->
registerParam
(
"conductivity_variation"
,
conductivity_variation
,
0.
,
_pat_parsmod
);
this
->
registerParam
(
"temperature_reference"
,
T_ref
,
0.
,
_pat_parsmod
);
this
->
registerParam
(
"capacity"
,
capacity
,
_pat_parsmod
);
this
->
registerParam
(
"density"
,
density
,
_pat_parsmod
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
initModel
()
{
getFEEngine
().
initShapeFunctions
(
_not_ghost
);
getFEEngine
().
initShapeFunctions
(
_ghost
);
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
initParallel
(
MeshPartition
*
partition
,
DataAccessor
*
data_accessor
)
{
AKANTU_DEBUG_IN
();
if
(
data_accessor
==
NULL
)
data_accessor
=
this
;
Synchronizer
&
synch_parallel
=
createParallelSynch
(
partition
,
data_accessor
);
synch_registry
->
registerSynchronizer
(
synch_parallel
,
_gst_htm_capacity
);
synch_registry
->
registerSynchronizer
(
synch_parallel
,
_gst_htm_temperature
);
synch_registry
->
registerSynchronizer
(
synch_parallel
,
_gst_htm_gradient_temperature
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
initPBC
()
{
AKANTU_DEBUG_IN
();
Model
::
initPBC
();
pbc_synch
=
new
PBCSynchronizer
(
pbc_pair
);
synch_registry
->
registerSynchronizer
(
*
pbc_synch
,
_gst_htm_capacity
);
synch_registry
->
registerSynchronizer
(
*
pbc_synch
,
_gst_htm_temperature
);
changeLocalEquationNumberForPBC
(
pbc_pair
,
1
);
// 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
();
while
(
it
!=
end
)
{
(
*
blocked_dofs
)((
*
it
).
first
,
0
)
=
true
;
++
it
;
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
initArrays
()
{
AKANTU_DEBUG_IN
();
UInt
nb_nodes
=
getFEEngine
().
getMesh
().
getNbNodes
();
std
::
stringstream
sstr_temp
;
sstr_temp
<<
Model
::
id
<<
":temperature"
;
std
::
stringstream
sstr_temp_rate
;
sstr_temp_rate
<<
Model
::
id
<<
":temperature_rate"
;
std
::
stringstream
sstr_inc
;
sstr_inc
<<
Model
::
id
<<
":increment"
;
std
::
stringstream
sstr_ext_flx
;
sstr_ext_flx
<<
Model
::
id
<<
":external_flux"
;
std
::
stringstream
sstr_residual
;
sstr_residual
<<
Model
::
id
<<
":residual"
;
std
::
stringstream
sstr_lump
;
sstr_lump
<<
Model
::
id
<<
":lumped"
;
std
::
stringstream
sstr_boun
;
sstr_boun
<<
Model
::
id
<<
":blocked_dofs"
;
temperature
=
&
(
alloc
<
Real
>
(
sstr_temp
.
str
(),
nb_nodes
,
1
,
REAL_INIT_VALUE
));
temperature_rate
=
&
(
alloc
<
Real
>
(
sstr_temp_rate
.
str
(),
nb_nodes
,
1
,
REAL_INIT_VALUE
));
increment
=
&
(
alloc
<
Real
>
(
sstr_inc
.
str
(),
nb_nodes
,
1
,
REAL_INIT_VALUE
));
external_heat_rate
=
&
(
alloc
<
Real
>
(
sstr_ext_flx
.
str
(),
nb_nodes
,
1
,
REAL_INIT_VALUE
));
residual
=
&
(
alloc
<
Real
>
(
sstr_residual
.
str
(),
nb_nodes
,
1
,
REAL_INIT_VALUE
));
capacity_lumped
=
&
(
alloc
<
Real
>
(
sstr_lump
.
str
(),
nb_nodes
,
1
,
REAL_INIT_VALUE
));
blocked_dofs
=
&
(
alloc
<
bool
>
(
sstr_boun
.
str
(),
nb_nodes
,
1
,
false
));
Mesh
::
ConnectivityTypeList
::
const_iterator
it
;
/* -------------------------------------------------------------------------- */
// byelementtype vectors
getFEEngine
().
getMesh
().
initElementTypeMapArray
(
temperature_on_qpoints
,
1
,
spatial_dimension
);
getFEEngine
().
getMesh
().
initElementTypeMapArray
(
temperature_gradient
,
spatial_dimension
,
spatial_dimension
);
getFEEngine
().
getMesh
().
initElementTypeMapArray
(
conductivity_on_qpoints
,
spatial_dimension
*
spatial_dimension
,
spatial_dimension
);
getFEEngine
().
getMesh
().
initElementTypeMapArray
(
k_gradt_on_qpoints
,
spatial_dimension
,
spatial_dimension
);
getFEEngine
().
getMesh
().
initElementTypeMapArray
(
bt_k_gT
,
1
,
spatial_dimension
,
true
);
getFEEngine
().
getMesh
().
initElementTypeMapArray
(
int_bt_k_gT
,
1
,
spatial_dimension
,
true
);
getFEEngine
().
getMesh
().
initElementTypeMapArray
(
thermal_energy
,
1
,
spatial_dimension
);
for
(
UInt
g
=
_not_ghost
;
g
<=
_ghost
;
++
g
)
{
GhostType
gt
=
(
GhostType
)
g
;
const
Mesh
::
ConnectivityTypeList
&
type_list
=
getFEEngine
().
getMesh
().
getConnectivityTypeList
(
gt
);
for
(
it
=
type_list
.
begin
();
it
!=
type_list
.
end
();
++
it
)
{
if
(
Mesh
::
getSpatialDimension
(
*
it
)
!=
spatial_dimension
)
continue
;
UInt
nb_element
=
getFEEngine
().
getMesh
().
getNbElement
(
*
it
,
gt
);
UInt
nb_quad_points
=
this
->
getFEEngine
().
getNbIntegrationPoints
(
*
it
,
gt
)
*
nb_element
;
temperature_on_qpoints
(
*
it
,
gt
).
resize
(
nb_quad_points
);
temperature_on_qpoints
(
*
it
,
gt
).
clear
();
temperature_gradient
(
*
it
,
gt
).
resize
(
nb_quad_points
);
temperature_gradient
(
*
it
,
gt
).
clear
();
conductivity_on_qpoints
(
*
it
,
gt
).
resize
(
nb_quad_points
);
conductivity_on_qpoints
(
*
it
,
gt
).
clear
();
k_gradt_on_qpoints
(
*
it
,
gt
).
resize
(
nb_quad_points
);
k_gradt_on_qpoints
(
*
it
,
gt
).
clear
();
bt_k_gT
(
*
it
,
gt
).
resize
(
nb_quad_points
);
bt_k_gT
(
*
it
,
gt
).
clear
();
int_bt_k_gT
(
*
it
,
gt
).
resize
(
nb_element
);
int_bt_k_gT
(
*
it
,
gt
).
clear
();
thermal_energy
(
*
it
,
gt
).
resize
(
nb_element
);
thermal_energy
(
*
it
,
gt
).
clear
();
}
}
/* -------------------------------------------------------------------------- */
dof_synchronizer
=
new
DOFSynchronizer
(
getFEEngine
().
getMesh
(),
1
);
dof_synchronizer
->
initLocalDOFEquationNumbers
();
dof_synchronizer
->
initGlobalDOFEquationNumbers
();
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
initSolver
(
__attribute__
((
unused
))
SolverOptions
&
options
)
{
#if !defined(AKANTU_USE_MUMPS)
// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake
AKANTU_DEBUG_ERROR
(
"You should at least activate one solver."
);
#else
UInt
nb_global_nodes
=
mesh
.
getNbGlobalNodes
();
delete
jacobian_matrix
;
std
::
stringstream
sstr
;
sstr
<<
Memory
::
id
<<
":jacobian_matrix"
;
jacobian_matrix
=
new
SparseMatrix
(
nb_global_nodes
,
_symmetric
,
sstr
.
str
(),
memory_id
);
jacobian_matrix
->
buildProfile
(
mesh
,
*
dof_synchronizer
,
1
);
delete
conductivity_matrix
;
std
::
stringstream
sstr_sti
;
sstr_sti
<<
Memory
::
id
<<
":conductivity_matrix"
;
conductivity_matrix
=
new
SparseMatrix
(
*
jacobian_matrix
,
sstr_sti
.
str
(),
memory_id
);
#ifdef AKANTU_USE_MUMPS
std
::
stringstream
sstr_solv
;
sstr_solv
<<
Memory
::
id
<<
":solver"
;
solver
=
new
SolverMumps
(
*
jacobian_matrix
,
sstr_solv
.
str
());
dof_synchronizer
->
initScatterGatherCommunicationScheme
();
#else
AKANTU_DEBUG_ERROR
(
"You should at least activate one solver."
);
#endif
//AKANTU_USE_MUMPS
if
(
solver
)
solver
->
initialize
(
options
);
#endif
//AKANTU_HAS_SOLVER
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
initImplicit
(
bool
dynamic
,
SolverOptions
&
solver_options
)
{
AKANTU_DEBUG_IN
();
method
=
dynamic
?
_implicit_dynamic
:
_static
;
initSolver
(
solver_options
);
if
(
method
==
_implicit_dynamic
)
{
if
(
integrator
)
delete
integrator
;
integrator
=
new
TrapezoidalRule1
();
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
HeatTransferModel
::~
HeatTransferModel
()
{
AKANTU_DEBUG_IN
();
if
(
integrator
)
delete
integrator
;
if
(
conductivity_matrix
)
delete
conductivity_matrix
;
if
(
capacity_matrix
)
delete
capacity_matrix
;
if
(
jacobian_matrix
)
delete
jacobian_matrix
;
if
(
solver
)
delete
solver
;
delete
pbc_synch
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
assembleCapacityLumped
(
const
GhostType
&
ghost_type
)
{
AKANTU_DEBUG_IN
();
FEEngine
&
fem
=
getFEEngine
();
const
Mesh
::
ConnectivityTypeList
&
type_list
=
fem
.
getMesh
().
getConnectivityTypeList
(
ghost_type
);
Mesh
::
ConnectivityTypeList
::
const_iterator
it
;
for
(
it
=
type_list
.
begin
();
it
!=
type_list
.
end
();
++
it
)
{
if
(
Mesh
::
getSpatialDimension
(
*
it
)
!=
spatial_dimension
)
continue
;
UInt
nb_element
=
getFEEngine
().
getMesh
().
getNbElement
(
*
it
,
ghost_type
);
UInt
nb_quadrature_points
=
getFEEngine
().
getNbIntegrationPoints
(
*
it
,
ghost_type
);
Array
<
Real
>
rho_1
(
nb_element
*
nb_quadrature_points
,
1
,
capacity
*
density
);
fem
.
assembleFieldLumped
(
rho_1
,
1
,
*
capacity_lumped
,
dof_synchronizer
->
getLocalDOFEquationNumbers
(),
*
it
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
assembleCapacityLumped
()
{
AKANTU_DEBUG_IN
();
capacity_lumped
->
clear
();
assembleCapacityLumped
(
_not_ghost
);
assembleCapacityLumped
(
_ghost
);
getSynchronizerRegistry
().
synchronize
(
_gst_htm_capacity
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
updateResidual
(
bool
compute_conductivity
)
{
AKANTU_DEBUG_IN
();
/// @f$ r = q_{ext} - q_{int} - C \dot T @f$
// start synchronization
synch_registry
->
asynchronousSynchronize
(
_gst_htm_temperature
);
// finalize communications
synch_registry
->
waitEndSynchronize
(
_gst_htm_temperature
);
//clear the array
/// first @f$ r = q_{ext} @f$
// residual->clear();
residual
->
copy
(
*
external_heat_rate
);
/// then @f$ r -= q_{int} @f$
// update the not ghost ones
updateResidual
(
_not_ghost
);
// update for the received ghosts
updateResidual
(
_ghost
);
/* if (method == _explicit_lumped_capacity) {
this->solveExplicitLumped();
}*/
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
assembleConductivityMatrix
(
bool
compute_conductivity
)
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_INFO
(
"Assemble the new stiffness matrix."
);
conductivity_matrix
->
clear
();
switch
(
mesh
.
getSpatialDimension
())
{
case
1
:
this
->
assembleConductivityMatrix
<
1
>
(
_not_ghost
,
compute_conductivity
);
break
;
case
2
:
this
->
assembleConductivityMatrix
<
2
>
(
_not_ghost
,
compute_conductivity
);
break
;
case
3
:
this
->
assembleConductivityMatrix
<
3
>
(
_not_ghost
,
compute_conductivity
);
break
;
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
HeatTransferModel
::
assembleConductivityMatrix
(
const
GhostType
&
ghost_type
,
bool
compute_conductivity
)
{
AKANTU_DEBUG_IN
();
Mesh
&
mesh
=
this
->
getFEEngine
().
getMesh
();
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
,
ghost_type
);
Mesh
::
type_iterator
last_type
=
mesh
.
lastType
(
spatial_dimension
,
ghost_type
);
for
(;
it
!=
last_type
;
++
it
)
{
this
->
assembleConductivityMatrix
<
dim
>
(
*
it
,
ghost_type
,
compute_conductivity
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template
<
UInt
dim
>
void
HeatTransferModel
::
assembleConductivityMatrix
(
const
ElementType
&
type
,
const
GhostType
&
ghost_type
,
bool
compute_conductivity
)
{
AKANTU_DEBUG_IN
();
SparseMatrix
&
K
=
*
conductivity_matrix
;
const
Array
<
Real
>
&
shapes_derivatives
=
this
->
getFEEngine
().
getShapesDerivatives
(
type
,
ghost_type
);
UInt
nb_element
=
mesh
.
getNbElement
(
type
,
ghost_type
);
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
type
);
UInt
nb_quadrature_points
=
getFEEngine
().
getNbIntegrationPoints
(
type
,
ghost_type
);
/// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
UInt
bt_d_b_size
=
nb_nodes_per_element
;
Array
<
Real
>
*
bt_d_b
=
new
Array
<
Real
>
(
nb_element
*
nb_quadrature_points
,
bt_d_b_size
*
bt_d_b_size
,
"B^t*D*B"
);
Matrix
<
Real
>
Bt_D
(
nb_nodes_per_element
,
dim
);
Array
<
Real
>::
const_iterator
<
Matrix
<
Real
>
>
shapes_derivatives_it
=
shapes_derivatives
.
begin
(
dim
,
nb_nodes_per_element
);
Array
<
Real
>::
iterator
<
Matrix
<
Real
>
>
Bt_D_B_it
=
bt_d_b
->
begin
(
bt_d_b_size
,
bt_d_b_size
);
if
(
compute_conductivity
)
this
->
computeConductivityOnQuadPoints
(
ghost_type
);
Array
<
Real
>::
iterator
<
Matrix
<
Real
>
>
D_it
=
conductivity_on_qpoints
(
type
,
ghost_type
).
begin
(
dim
,
dim
);
Array
<
Real
>::
iterator
<
Matrix
<
Real
>
>
D_end
=
conductivity_on_qpoints
(
type
,
ghost_type
).
end
(
dim
,
dim
);
for
(;
D_it
!=
D_end
;
++
D_it
,
++
Bt_D_B_it
,
++
shapes_derivatives_it
)
{
Matrix
<
Real
>
&
D
=
*
D_it
;
const
Matrix
<
Real
>
&
B
=
*
shapes_derivatives_it
;
Matrix
<
Real
>
&
Bt_D_B
=
*
Bt_D_B_it
;
Bt_D
.
mul
<
true
,
false
>
(
B
,
D
);
Bt_D_B
.
mul
<
false
,
false
>
(
Bt_D
,
B
);
}
/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
Array
<
Real
>
*
K_e
=
new
Array
<
Real
>
(
nb_element
,
bt_d_b_size
*
bt_d_b_size
,
"K_e"
);
this
->
getFEEngine
().
integrate
(
*
bt_d_b
,
*
K_e
,
bt_d_b_size
*
bt_d_b_size
,
type
,
ghost_type
);
delete
bt_d_b
;
this
->
getFEEngine
().
assembleMatrix
(
*
K_e
,
K
,
1
,
type
,
ghost_type
);
delete
K_e
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
updateResidualInternal
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_INFO
(
"Update the residual"
);
// f = q_ext - q_int - Mddot = q - Mddot;
if
(
method
!=
_static
)
{
// f -= Mddot
if
(
capacity_matrix
)
{
// if full mass_matrix
Array
<
Real
>
*
Mddot
=
new
Array
<
Real
>
(
*
temperature_rate
,
true
,
"Mddot"
);
*
Mddot
*=
*
capacity_matrix
;
*
residual
-=
*
Mddot
;
delete
Mddot
;
}
else
if
(
capacity_lumped
)
{
// else lumped mass
UInt
nb_nodes
=
temperature_rate
->
getSize
();
UInt
nb_degree_of_freedom
=
temperature_rate
->
getNbComponent
();
Real
*
capacity_val
=
capacity_lumped
->
storage
();
Real
*
temp_rate_val
=
temperature_rate
->
storage
();
Real
*
res_val
=
residual
->
storage
();
bool
*
blocked_dofs_val
=
blocked_dofs
->
storage
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
*
nb_degree_of_freedom
;
++
n
)
{
if
(
!
(
*
blocked_dofs_val
))
{
*
res_val
-=
*
temp_rate_val
*
*
capacity_val
;
}
blocked_dofs_val
++
;
res_val
++
;
capacity_val
++
;
temp_rate_val
++
;
}
}
else
{
AKANTU_DEBUG_ERROR
(
"No function called to assemble the mass matrix."
);
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
solveStatic
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_INFO
(
"Solving Ku = f"
);
AKANTU_DEBUG_ASSERT
(
conductivity_matrix
!=
NULL
,
"You should first initialize the implicit solver and assemble the stiffness matrix"
);
UInt
nb_nodes
=
temperature
->
getSize
();
UInt
nb_degree_of_freedom
=
temperature
->
getNbComponent
()
*
nb_nodes
;
jacobian_matrix
->
copyContent
(
*
conductivity_matrix
);
jacobian_matrix
->
applyBoundary
(
*
blocked_dofs
);
increment
->
clear
();
solver
->
setRHS
(
*
residual
);
solver
->
factorize
();
solver
->
solve
(
*
increment
);
Real
*
increment_val
=
increment
->
storage
();
Real
*
temperature_val
=
temperature
->
storage
();
bool
*
blocked_dofs_val
=
blocked_dofs
->
storage
();
for
(
UInt
j
=
0
;
j
<
nb_degree_of_freedom
;
++
j
,
++
temperature_val
,
++
increment_val
,
++
blocked_dofs_val
)
{
if
(
!
(
*
blocked_dofs_val
))
{
*
temperature_val
+=
*
increment_val
;
}
else
{
*
increment_val
=
0.0
;
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
computeConductivityOnQuadPoints
(
const
GhostType
&
ghost_type
)
{
const
Mesh
::
ConnectivityTypeList
&
type_list
=
this
->
getFEEngine
().
getMesh
().
getConnectivityTypeList
(
ghost_type
);
Mesh
::
ConnectivityTypeList
::
const_iterator
it
;
for
(
it
=
type_list
.
begin
();
it
!=
type_list
.
end
();
++
it
)
{
if
(
Mesh
::
getSpatialDimension
(
*
it
)
!=
spatial_dimension
)
continue
;
Array
<
Real
>
&
temperature_interpolated
=
temperature_on_qpoints
(
*
it
,
ghost_type
);
//compute the temperature on quadrature points
this
->
getFEEngine
().
interpolateOnIntegrationPoints
(
*
temperature
,
temperature_interpolated
,
1
,
*
it
,
ghost_type
);
Array
<
Real
>::
matrix_iterator
C_it
=
conductivity_on_qpoints
(
*
it
,
ghost_type
).
begin
(
spatial_dimension
,
spatial_dimension
);
Array
<
Real
>::
matrix_iterator
C_end
=
conductivity_on_qpoints
(
*
it
,
ghost_type
).
end
(
spatial_dimension
,
spatial_dimension
);
Array
<
Real
>::
iterator
<
Real
>
T_it
=
temperature_interpolated
.
begin
();
for
(;
C_it
!=
C_end
;
++
C_it
,
++
T_it
)
{
Matrix
<
Real
>
&
C
=
*
C_it
;
Real
&
T
=
*
T_it
;
C
=
conductivity
;
Matrix
<
Real
>
variation
(
spatial_dimension
,
spatial_dimension
,
conductivity_variation
*
(
T
-
T_ref
));
C
+=
conductivity_variation
;
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
computeKgradT
(
const
GhostType
&
ghost_type
,
bool
compute_conductivity
)
{
if
(
compute_conductivity
)
computeConductivityOnQuadPoints
(
ghost_type
);
const
Mesh
::
ConnectivityTypeList
&
type_list
=
this
->
getFEEngine
().
getMesh
().
getConnectivityTypeList
(
ghost_type
);
Mesh
::
ConnectivityTypeList
::
const_iterator
it
;
for
(
it
=
type_list
.
begin
();
it
!=
type_list
.
end
();
++
it
)
{
const
ElementType
&
type
=
*
it
;
if
(
Mesh
::
getSpatialDimension
(
*
it
)
!=
spatial_dimension
)
continue
;
Array
<
Real
>
&
gradient
=
temperature_gradient
(
*
it
,
ghost_type
);
this
->
getFEEngine
().
gradientOnIntegrationPoints
(
*
temperature
,
gradient
,
1
,
*
it
,
ghost_type
);
Array
<
Real
>::
matrix_iterator
C_it
=
conductivity_on_qpoints
(
*
it
,
ghost_type
).
begin
(
spatial_dimension
,
spatial_dimension
);
Array
<
Real
>::
vector_iterator
BT_it
=
gradient
.
begin
(
spatial_dimension
);
Array
<
Real
>::
vector_iterator
k_BT_it
=
k_gradt_on_qpoints
(
type
,
ghost_type
).
begin
(
spatial_dimension
);
Array
<
Real
>::
vector_iterator
k_BT_end
=
k_gradt_on_qpoints
(
type
,
ghost_type
).
end
(
spatial_dimension
);
for
(;
k_BT_it
!=
k_BT_end
;
++
k_BT_it
,
++
BT_it
,
++
C_it
)
{
Vector
<
Real
>
&
k_BT
=
*
k_BT_it
;
Vector
<
Real
>
&
BT
=
*
BT_it
;
Matrix
<
Real
>
&
C
=
*
C_it
;
k_BT
.
mul
<
false
>
(
C
,
BT
);
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
updateResidual
(
const
GhostType
&
ghost_type
,
bool
compute_conductivity
)
{
AKANTU_DEBUG_IN
();
const
Mesh
::
ConnectivityTypeList
&
type_list
=
this
->
getFEEngine
().
getMesh
().
getConnectivityTypeList
(
ghost_type
);
Mesh
::
ConnectivityTypeList
::
const_iterator
it
;
for
(
it
=
type_list
.
begin
();
it
!=
type_list
.
end
();
++
it
)
{
if
(
Mesh
::
getSpatialDimension
(
*
it
)
!=
spatial_dimension
)
continue
;
Array
<
Real
>
&
shapes_derivatives
=
const_cast
<
Array
<
Real
>
&>
(
getFEEngine
().
getShapesDerivatives
(
*
it
,
ghost_type
));
UInt
nb_nodes_per_element
=
Mesh
::
getNbNodesPerElement
(
*
it
);
// compute k \grad T
computeKgradT
(
ghost_type
,
compute_conductivity
);
Array
<
Real
>::
vector_iterator
k_BT_it
=
k_gradt_on_qpoints
(
*
it
,
ghost_type
).
begin
(
spatial_dimension
);
Array
<
Real
>::
matrix_iterator
B_it
=
shapes_derivatives
.
begin
(
spatial_dimension
,
nb_nodes_per_element
);
Array
<
Real
>::
vector_iterator
Bt_k_BT_it
=
bt_k_gT
(
*
it
,
ghost_type
).
begin
(
nb_nodes_per_element
);
Array
<
Real
>::
vector_iterator
Bt_k_BT_end
=
bt_k_gT
(
*
it
,
ghost_type
).
end
(
nb_nodes_per_element
);
for
(;
Bt_k_BT_it
!=
Bt_k_BT_end
;
++
Bt_k_BT_it
,
++
B_it
,
++
k_BT_it
)
{
Vector
<
Real
>
&
k_BT
=
*
k_BT_it
;
Vector
<
Real
>
&
Bt_k_BT
=
*
Bt_k_BT_it
;
Matrix
<
Real
>
&
B
=
*
B_it
;
Bt_k_BT
.
mul
<
true
>
(
B
,
k_BT
);
}
this
->
getFEEngine
().
integrate
(
bt_k_gT
(
*
it
,
ghost_type
),
int_bt_k_gT
(
*
it
,
ghost_type
),
nb_nodes_per_element
,
*
it
,
ghost_type
);
this
->
getFEEngine
().
assembleArray
(
int_bt_k_gT
(
*
it
,
ghost_type
),
*
residual
,
dof_synchronizer
->
getLocalDOFEquationNumbers
(),
1
,
*
it
,
ghost_type
,
empty_filter
,
-
1
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
solveExplicitLumped
()
{
AKANTU_DEBUG_IN
();
/// finally @f$ r -= C \dot T @f$
// lumped C
UInt
nb_nodes
=
temperature_rate
->
getSize
();
UInt
nb_degree_of_freedom
=
temperature_rate
->
getNbComponent
();
Real
*
capacity_val
=
capacity_lumped
->
storage
();
Real
*
temp_rate_val
=
temperature_rate
->
storage
();
Real
*
res_val
=
residual
->
storage
();
bool
*
blocked_dofs_val
=
blocked_dofs
->
storage
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
*
nb_degree_of_freedom
;
++
n
)
{
if
(
!
(
*
blocked_dofs_val
))
{
*
res_val
-=
*
capacity_val
*
*
temp_rate_val
;
}
blocked_dofs_val
++
;
res_val
++
;
capacity_val
++
;
temp_rate_val
++
;
}
#ifndef AKANTU_NDEBUG
getSynchronizerRegistry
().
synchronize
(
akantu
::
_gst_htm_gradient_temperature
);
#endif
capacity_val
=
capacity_lumped
->
storage
();
res_val
=
residual
->
storage
();
blocked_dofs_val
=
blocked_dofs
->
storage
();
Real
*
inc
=
increment
->
storage
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
*
nb_degree_of_freedom
;
++
n
)
{
if
(
!
(
*
blocked_dofs_val
))
{
*
inc
=
(
*
res_val
/
*
capacity_val
);
}
res_val
++
;
blocked_dofs_val
++
;
inc
++
;
capacity_val
++
;
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
explicitPred
()
{
AKANTU_DEBUG_IN
();
AKANTU_DEBUG_ASSERT
(
integrator
,
"integrator was not instanciated"
);
integrator
->
integrationSchemePred
(
time_step
,
*
temperature
,
*
temperature_rate
,
*
blocked_dofs
);
UInt
nb_nodes
=
temperature
->
getSize
();
UInt
nb_degree_of_freedom
=
temperature
->
getNbComponent
();
Real
*
temp
=
temperature
->
storage
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
*
nb_degree_of_freedom
;
++
n
,
++
temp
)
if
(
*
temp
<
0.
)
*
temp
=
0.
;
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
explicitCorr
()
{
AKANTU_DEBUG_IN
();
integrator
->
integrationSchemeCorrTempRate
(
time_step
,
*
temperature
,
*
temperature_rate
,
*
blocked_dofs
,
*
increment
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
implicitPred
(){
AKANTU_DEBUG_IN
();
if
(
method
==
_implicit_dynamic
)
integrator
->
integrationSchemePred
(
time_step
,
*
temperature
,
*
temperature_rate
,
*
blocked_dofs
);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
implicitCorr
(){
AKANTU_DEBUG_IN
();
if
(
method
==
_implicit_dynamic
)
{
integrator
->
integrationSchemeCorrTemp
(
time_step
,
*
temperature
,
*
temperature_rate
,
*
blocked_dofs
,
*
increment
);
}
else
{
UInt
nb_nodes
=
temperature
->
getSize
();
UInt
nb_degree_of_freedom
=
temperature
->
getNbComponent
()
*
nb_nodes
;
Real
*
incr_val
=
increment
->
storage
();
Real
*
temp_val
=
temperature
->
storage
();
bool
*
boun_val
=
blocked_dofs
->
storage
();
for
(
UInt
j
=
0
;
j
<
nb_degree_of_freedom
;
++
j
,
++
temp_val
,
++
incr_val
,
++
boun_val
){
*
incr_val
*=
(
1.
-
*
boun_val
);
*
temp_val
+=
*
incr_val
;
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
Real
HeatTransferModel
::
getStableTimeStep
()
{
AKANTU_DEBUG_IN
();
Real
el_size
;
Real
min_el_size
=
std
::
numeric_limits
<
Real
>::
max
();
Real
conductivitymax
=
conductivity
(
0
,
0
);
//get the biggest parameter from k11 until k33//
for
(
UInt
i
=
0
;
i
<
spatial_dimension
;
i
++
)
for
(
UInt
j
=
0
;
j
<
spatial_dimension
;
j
++
)
conductivitymax
=
std
::
max
(
conductivity
(
i
,
j
),
conductivitymax
);
const
Mesh
::
ConnectivityTypeList
&
type_list
=
getFEEngine
().
getMesh
().
getConnectivityTypeList
();
Mesh
::
ConnectivityTypeList
::
const_iterator
it
;
for
(
it
=
type_list
.
begin
();
it
!=
type_list
.
end
();
++
it
)
{
if
(
getFEEngine
().
getMesh
().
getSpatialDimension
(
*
it
)
!=
spatial_dimension
)
continue
;
UInt
nb_nodes_per_element
=
getFEEngine
().
getMesh
().
getNbNodesPerElement
(
*
it
);
Array
<
Real
>
coord
(
0
,
nb_nodes_per_element
*
spatial_dimension
);
FEEngine
::
extractNodalToElementField
(
getFEEngine
().
getMesh
(),
getFEEngine
().
getMesh
().
getNodes
(),
coord
,
*
it
,
_not_ghost
);
Array
<
Real
>::
matrix_iterator
el_coord
=
coord
.
begin
(
spatial_dimension
,
nb_nodes_per_element
);
UInt
nb_element
=
getFEEngine
().
getMesh
().
getNbElement
(
*
it
);
for
(
UInt
el
=
0
;
el
<
nb_element
;
++
el
,
++
el_coord
)
{
el_size
=
getFEEngine
().
getElementInradius
(
*
el_coord
,
*
it
);
min_el_size
=
std
::
min
(
min_el_size
,
el_size
);
}
AKANTU_DEBUG_INFO
(
"The minimum element size : "
<<
min_el_size
<<
" and the max conductivity is : "
<<
conductivitymax
);
}
Real
min_dt
=
2
*
min_el_size
*
min_el_size
*
density
*
capacity
/
conductivitymax
;
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
min_dt
,
1
,
_so_min
);
AKANTU_DEBUG_OUT
();
return
min_dt
;
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
readMaterials
()
{
std
::
pair
<
Parser
::
const_section_iterator
,
Parser
::
const_section_iterator
>
sub_sect
=
this
->
parser
->
getSubSections
(
_st_heat
);
Parser
::
const_section_iterator
it
=
sub_sect
.
first
;
const
ParserSection
&
section
=
*
it
;
this
->
parseSection
(
section
);
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
initFull
(
const
ModelOptions
&
options
){
Model
::
initFull
(
options
);
readMaterials
();
const
HeatTransferModelOptions
&
my_options
=
dynamic_cast
<
const
HeatTransferModelOptions
&>
(
options
);
method
=
my_options
.
analysis_method
;
//initialize the vectors
initArrays
();
temperature
->
clear
();
temperature_rate
->
clear
();
external_heat_rate
->
clear
();
// initialize pbc
if
(
pbc_pair
.
size
()
!=
0
)
initPBC
();
if
(
method
==
_explicit_lumped_capacity
){
integrator
=
new
ForwardEuler
();
}
if
(
method
==
_implicit_dynamic
)
{
initImplicit
(
true
);
}
if
(
method
==
_static
)
{
initImplicit
(
false
);
}
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
assembleCapacity
()
{
AKANTU_DEBUG_IN
();
if
(
!
capacity_matrix
)
{
std
::
stringstream
sstr
;
sstr
<<
id
<<
":capacity_matrix"
;
capacity_matrix
=
new
SparseMatrix
(
*
jacobian_matrix
,
sstr
.
str
(),
memory_id
);
}
assembleCapacity
(
_not_ghost
);
// assembleMass(_ghost);
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
assembleCapacity
(
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
MyFEEngineType
&
fem
=
getFEEngineClass
<
MyFEEngineType
>
();
Array
<
Real
>
rho_1
(
0
,
1
);
//UInt nb_element;
capacity_matrix
->
clear
();
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
;
computeRho
(
rho_1
,
type
,
ghost_type
);
fem
.
assembleFieldMatrix
(
rho_1
,
1
,
*
capacity_matrix
,
type
,
ghost_type
);
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
computeRho
(
Array
<
Real
>
&
rho
,
ElementType
type
,
GhostType
ghost_type
)
{
AKANTU_DEBUG_IN
();
FEEngine
&
fem
=
this
->
getFEEngine
();
UInt
nb_element
=
fem
.
getMesh
().
getNbElement
(
type
,
ghost_type
);
UInt
nb_quadrature_points
=
fem
.
getNbIntegrationPoints
(
type
,
ghost_type
);
rho
.
resize
(
nb_element
*
nb_quadrature_points
);
Real
*
rho_1_val
=
rho
.
storage
();
/// compute @f$ rho @f$ for each nodes of each element
for
(
UInt
el
=
0
;
el
<
nb_element
;
++
el
)
{
for
(
UInt
n
=
0
;
n
<
nb_quadrature_points
;
++
n
)
{
*
rho_1_val
++
=
this
->
capacity
;
}
}
AKANTU_DEBUG_OUT
();
}
/* -------------------------------------------------------------------------- */
template
<>
bool
HeatTransferModel
::
testConvergence
<
_scc_increment
>
(
Real
tolerance
,
Real
&
error
){
AKANTU_DEBUG_IN
();
UInt
nb_nodes
=
temperature
->
getSize
();
UInt
nb_degree_of_freedom
=
temperature
->
getNbComponent
();
error
=
0
;
Real
norm
[
2
]
=
{
0.
,
0.
};
Real
*
increment_val
=
increment
->
storage
();
bool
*
blocked_dofs_val
=
blocked_dofs
->
storage
();
Real
*
temperature_val
=
temperature
->
storage
();
for
(
UInt
n
=
0
;
n
<
nb_nodes
;
++
n
)
{
bool
is_local_node
=
mesh
.
isLocalOrMasterNode
(
n
);
for
(
UInt
d
=
0
;
d
<
nb_degree_of_freedom
;
++
d
)
{
if
(
!
(
*
blocked_dofs_val
)
&&
is_local_node
)
{
norm
[
0
]
+=
*
increment_val
*
*
increment_val
;
norm
[
1
]
+=
*
temperature_val
*
*
temperature_val
;
}
blocked_dofs_val
++
;
increment_val
++
;
temperature_val
++
;
}
}
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
norm
,
2
,
_so_sum
);
norm
[
0
]
=
sqrt
(
norm
[
0
]);
norm
[
1
]
=
sqrt
(
norm
[
1
]);
AKANTU_DEBUG_ASSERT
(
!
Math
::
isnan
(
norm
[
0
]),
"Something goes wrong in the solve phase"
);
if
(
norm
[
1
]
<
Math
::
getTolerance
())
{
error
=
norm
[
0
];
AKANTU_DEBUG_OUT
();
// cout<<"Error 1: "<<error<<endl;
return
error
<
tolerance
;
}
AKANTU_DEBUG_OUT
();
if
(
norm
[
1
]
>
Math
::
getTolerance
())
error
=
norm
[
0
]
/
norm
[
1
];
else
error
=
norm
[
0
];
//In case the total displacement is zero!
// cout<<"Error 2: "<<error<<endl;
return
(
error
<
tolerance
);
}
/* -------------------------------------------------------------------------- */
void
HeatTransferModel
::
initFEEngineBoundary
(
bool
create_surface
)
{
if
(
create_surface
)
MeshUtils
::
buildFacets
(
getFEEngine
().
getMesh
());
FEEngine
&
fem_boundary
=
getFEEngineBoundary
();
fem_boundary
.
initShapeFunctions
();
fem_boundary
.
computeNormalsOnIntegrationPoints
();
}
/* -------------------------------------------------------------------------- */
Real
HeatTransferModel
::
computeThermalEnergyByNode
()
{
AKANTU_DEBUG_IN
();
Real
ethermal
=
0.
;
Array
<
Real
>::
vector_iterator
heat_rate_it
=
residual
->
begin
(
residual
->
getNbComponent
());
Array
<
Real
>::
vector_iterator
heat_rate_end
=
residual
->
end
(
residual
->
getNbComponent
());
UInt
n
=
0
;
for
(;
heat_rate_it
!=
heat_rate_end
;
++
heat_rate_it
,
++
n
)
{
Real
heat
=
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
;
Vector
<
Real
>
&
heat_rate
=
*
heat_rate_it
;
for
(
UInt
i
=
0
;
i
<
heat_rate
.
size
();
++
i
)
{
if
(
count_node
)
heat
+=
heat_rate
[
i
]
*
time_step
;
}
ethermal
+=
heat
;
}
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
ethermal
,
1
,
_so_sum
);
AKANTU_DEBUG_OUT
();
return
ethermal
;
}
/* -------------------------------------------------------------------------- */
template
<
class
iterator
>
void
HeatTransferModel
::
getThermalEnergy
(
iterator
Eth
,
Array
<
Real
>::
const_iterator
<
Real
>
T_it
,
Array
<
Real
>::
const_iterator
<
Real
>
T_end
)
const
{
for
(;
T_it
!=
T_end
;
++
T_it
,
++
Eth
)
{
*
Eth
=
capacity
*
density
*
*
T_it
;
}
}
/* -------------------------------------------------------------------------- */
Real
HeatTransferModel
::
getThermalEnergy
(
const
ElementType
&
type
,
UInt
index
)
{
AKANTU_DEBUG_IN
();
UInt
nb_quadrature_points
=
getFEEngine
().
getNbIntegrationPoints
(
type
);
Vector
<
Real
>
Eth_on_quarature_points
(
nb_quadrature_points
);
Array
<
Real
>::
iterator
<
Real
>
T_it
=
this
->
temperature_on_qpoints
(
type
).
begin
();
T_it
+=
index
*
nb_quadrature_points
;
Array
<
Real
>::
iterator
<
Real
>
T_end
=
T_it
+
nb_quadrature_points
;
getThermalEnergy
(
Eth_on_quarature_points
.
storage
(),
T_it
,
T_end
);
return
getFEEngine
().
integrate
(
Eth_on_quarature_points
,
type
,
index
);
}
/* -------------------------------------------------------------------------- */
Real
HeatTransferModel
::
getThermalEnergy
()
{
Real
Eth
=
0
;
Mesh
&
mesh
=
getFEEngine
().
getMesh
();
Mesh
::
type_iterator
it
=
mesh
.
firstType
(
spatial_dimension
);
Mesh
::
type_iterator
last_type
=
mesh
.
lastType
(
spatial_dimension
);
for
(;
it
!=
last_type
;
++
it
)
{
UInt
nb_element
=
getFEEngine
().
getMesh
().
getNbElement
(
*
it
,
_not_ghost
);
UInt
nb_quadrature_points
=
getFEEngine
().
getNbIntegrationPoints
(
*
it
,
_not_ghost
);
Array
<
Real
>
Eth_per_quad
(
nb_element
*
nb_quadrature_points
,
1
);
Array
<
Real
>::
iterator
<
Real
>
T_it
=
this
->
temperature_on_qpoints
(
*
it
).
begin
();
Array
<
Real
>::
iterator
<
Real
>
T_end
=
this
->
temperature_on_qpoints
(
*
it
).
end
();
getThermalEnergy
(
Eth_per_quad
.
begin
(),
T_it
,
T_end
);
Eth
+=
getFEEngine
().
integrate
(
Eth_per_quad
,
*
it
);
}
return
Eth
;
}
/* -------------------------------------------------------------------------- */
Real
HeatTransferModel
::
getEnergy
(
const
std
::
string
&
id
)
{
AKANTU_DEBUG_IN
();
Real
energy
=
0
;
if
(
"thermal"
)
energy
=
getThermalEnergy
();
// reduction sum over all processors
StaticCommunicator
::
getStaticCommunicator
().
allReduce
(
&
energy
,
1
,
_so_sum
);
AKANTU_DEBUG_OUT
();
return
energy
;
}
/* -------------------------------------------------------------------------- */
Real
HeatTransferModel
::
getEnergy
(
const
std
::
string
&
energy_id
,
const
ElementType
&
type
,
UInt
index
)
{
AKANTU_DEBUG_IN
();
Real
energy
=
0.
;
if
(
"thermal"
)
energy
=
getThermalEnergy
(
type
,
index
);
AKANTU_DEBUG_OUT
();
return
energy
;
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
dumper
::
Field
*
HeatTransferModel
::
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
;
}
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
HeatTransferModel
::
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
[
"temperature"
]
=
temperature
;
real_nodal_fields
[
"temperature_rate"
]
=
temperature_rate
;
real_nodal_fields
[
"external_heat_rate"
]
=
external_heat_rate
;
real_nodal_fields
[
"residual"
]
=
residual
;
real_nodal_fields
[
"capacity_lumped"
]
=
capacity_lumped
;
real_nodal_fields
[
"increment"
]
=
increment
;
dumper
::
Field
*
field
=
mesh
.
createNodalField
(
real_nodal_fields
[
field_name
],
group_name
);
return
field
;
}
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
HeatTransferModel
::
createElementalField
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
,
const
UInt
&
spatial_dimension
,
const
ElementKind
&
element_kind
){
dumper
::
Field
*
field
=
NULL
;
if
(
field_name
==
"partitions"
)
field
=
mesh
.
createElementalField
<
UInt
,
dumper
::
ElementPartitionField
>
(
mesh
.
getConnectivities
(),
group_name
,
this
->
spatial_dimension
,
element_kind
);
else
if
(
field_name
==
"temperature_gradient"
){
ElementTypeMap
<
UInt
>
nb_data_per_elem
=
this
->
mesh
.
getNbDataPerElem
(
temperature_gradient
,
element_kind
);
field
=
mesh
.
createElementalField
<
Real
,
dumper
::
InternalMaterialField
>
(
temperature_gradient
,
group_name
,
this
->
spatial_dimension
,
element_kind
,
nb_data_per_elem
);
}
else
if
(
field_name
==
"conductivity"
){
ElementTypeMap
<
UInt
>
nb_data_per_elem
=
this
->
mesh
.
getNbDataPerElem
(
conductivity_on_qpoints
,
element_kind
);
field
=
mesh
.
createElementalField
<
Real
,
dumper
::
InternalMaterialField
>
(
conductivity_on_qpoints
,
group_name
,
this
->
spatial_dimension
,
element_kind
,
nb_data_per_elem
);
}
return
field
;
}
/* -------------------------------------------------------------------------- */
#else
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
HeatTransferModel
::
createElementalField
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
,
const
ElementKind
&
element_kind
){
return
NULL
;
}
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
HeatTransferModel
::
createNodalFieldBool
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
)
{
return
NULL
;
}
/* -------------------------------------------------------------------------- */
dumper
::
Field
*
HeatTransferModel
::
createNodalFieldReal
(
const
std
::
string
&
field_name
,
const
std
::
string
&
group_name
,
bool
padding_flag
)
{
return
NULL
;
}
#endif
__END_AKANTU__
Event Timeline
Log In to Comment