diff --git a/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index 548b671ac..e2e2c114e 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,1282 +1,1305 @@ /** * @file heat_transfer_model.cc * * @author Guillaume Anciaux * @author Lucas Frerot * @author David Simon Kammer * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @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 . * */ /* -------------------------------------------------------------------------- */ -#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 "aka_common.hh" +#include "aka_math.hh" +#include "dumpable_inline_impl.hh" #include "fe_engine_template.hh" +#include "generalized_trapezoidal.hh" +#include "group_manager_inline_impl.cc" #include "mesh.hh" -#include "static_communicator.hh" #include "parser.hh" -#include "generalized_trapezoidal.hh" +#include "static_communicator.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_elemental_field.hh" #include "dumper_internal_material_field.hh" +#include "dumper_paraview.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(sstr.str(), mesh, spatial_dimension); this->temperature = NULL; this->residual = NULL; this->blocked_dofs = NULL; #ifdef AKANTU_USE_IOHELPER - this->mesh.registerDumper("paraview_all", id, true); + this->mesh.registerDumper("heat_transfer", 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::iterator it = pbc_pair.begin(); std::map::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(sstr_temp.str(), nb_nodes, 1, REAL_INIT_VALUE)); temperature_rate = &(alloc(sstr_temp_rate.str(), nb_nodes, 1, REAL_INIT_VALUE)); increment = &(alloc(sstr_inc.str(), nb_nodes, 1, REAL_INIT_VALUE)); external_heat_rate = &(alloc(sstr_ext_flx.str(), nb_nodes, 1, REAL_INIT_VALUE)); residual = &(alloc(sstr_residual.str(), nb_nodes, 1, REAL_INIT_VALUE)); capacity_lumped = &(alloc(sstr_lump.str(), nb_nodes, 1, REAL_INIT_VALUE)); blocked_dofs = &(alloc(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 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( __attribute__((unused)) 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 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(*it, ghost_type, compute_conductivity); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template void HeatTransferModel::assembleConductivityMatrix(const ElementType & type, const GhostType & ghost_type, bool compute_conductivity) { AKANTU_DEBUG_IN(); SparseMatrix & K = *conductivity_matrix; const Array & 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 * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix Bt_D(nb_nodes_per_element, dim); Array::const_iterator > shapes_derivatives_it = shapes_derivatives.begin(dim, nb_nodes_per_element); Array::iterator > Bt_D_B_it = bt_d_b->begin(bt_d_b_size, bt_d_b_size); if (compute_conductivity) this->computeConductivityOnQuadPoints(ghost_type); Array::iterator > D_it = conductivity_on_qpoints(type, ghost_type).begin(dim, dim); Array::iterator > 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 & D = *D_it; const Matrix & B = *shapes_derivatives_it; Matrix & Bt_D_B = *Bt_D_B_it; Bt_D.mul(B, D); Bt_D_B.mul(Bt_D, B); } /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(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 * Mddot = new Array(*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; + *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 & 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::matrix_iterator C_it = conductivity_on_qpoints(*it, ghost_type) .begin(spatial_dimension, spatial_dimension); Array::matrix_iterator C_end = conductivity_on_qpoints(*it, ghost_type) .end(spatial_dimension, spatial_dimension); Array::iterator T_it = temperature_interpolated.begin(); for (; C_it != C_end; ++C_it, ++T_it) { Matrix & C = *C_it; Real & T = *T_it; C = conductivity; Matrix 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 & gradient = temperature_gradient(*it, ghost_type); this->getFEEngine().gradientOnIntegrationPoints(*temperature, gradient, 1, *it, ghost_type); Array::matrix_iterator C_it = conductivity_on_qpoints(*it, ghost_type) .begin(spatial_dimension, spatial_dimension); Array::vector_iterator BT_it = gradient.begin(spatial_dimension); Array::vector_iterator k_BT_it = k_gradt_on_qpoints(type, ghost_type).begin(spatial_dimension); Array::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 & k_BT = *k_BT_it; Vector & BT = *BT_it; Matrix & C = *C_it; k_BT.mul(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 & shapes_derivatives = const_cast &>( getFEEngine().getShapesDerivatives(*it, ghost_type)); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); // compute k \grad T computeKgradT(ghost_type, compute_conductivity); Array::vector_iterator k_BT_it = k_gradt_on_qpoints(*it, ghost_type).begin(spatial_dimension); Array::matrix_iterator B_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::vector_iterator Bt_k_BT_it = bt_k_gT(*it, ghost_type).begin(nb_nodes_per_element); Array::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 & k_BT = *k_BT_it; Vector & Bt_k_BT = *Bt_k_BT_it; Matrix & B = *B_it; Bt_k_BT.mul(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; + *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::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 coord(0, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(getFEEngine().getMesh(), getFEEngine().getMesh().getNodes(), coord, *it, _not_ghost); Array::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 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(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(); } /* -------------------------------------------------------------------------- */ class ComputeRhoFunctor { public: ComputeRhoFunctor(const HeatTransferModel & model) : model(model){}; void operator()(Matrix & rho, const Element &, const Matrix &) const { rho.set(model.getCapacity()); } private: const HeatTransferModel & model; }; /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacity(GhostType ghost_type) { AKANTU_DEBUG_IN(); MyFEEngineType & fem = getFEEngineClass(); ComputeRhoFunctor rho_functor(*this); - 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; fem.assembleFieldMatrix(rho_functor, 1, *capacity_matrix, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeRho(Array & 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: "< Math::getTolerance()) error = norm[0] / norm[1]; else error = norm[0]; // In case the total displacement is zero! // cout<<"Error 2: "<::vector_iterator heat_rate_it = residual->begin(residual->getNbComponent()); Array::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 & 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(ðermal, 1, _so_sum); AKANTU_DEBUG_OUT(); return ethermal; } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::getThermalEnergy( iterator Eth, Array::const_iterator T_it, Array::const_iterator T_end) const { for (; T_it != T_end; ++T_it, ++Eth) { - *Eth = capacity * density ** T_it; + *Eth = capacity * density * *T_it; } } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Vector Eth_on_quarature_points(nb_quadrature_points); Array::iterator T_it = this->temperature_on_qpoints(type).begin(); T_it += index * nb_quadrature_points; Array::iterator 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 Eth_per_quad(nb_element * nb_quadrature_points, 1); Array::iterator T_it = this->temperature_on_qpoints(*it).begin(); Array::iterator 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 (id == "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 & id, const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real energy = 0.; if (id == "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, __attribute__((unused)) bool padding_flag) { std::map *> 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, __attribute__((unused)) bool padding_flag) { std::map *> 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, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const UInt & spatial_dimension, const ElementKind & element_kind) { dumper::Field * field = NULL; if (field_name == "partitions") field = mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); else if (field_name == "temperature_gradient") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(temperature_gradient, element_kind); field = mesh.createElementalField( temperature_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "conductivity") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(conductivity_on_qpoints, element_kind); field = mesh.createElementalField( conductivity_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ dumper::Field * HeatTransferModel::createElementalField( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const ElementKind & element_kind) { return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * HeatTransferModel::createNodalFieldBool( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * HeatTransferModel::createNodalFieldReal( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return NULL; } #endif +void HeatTransferModel::dump(const std::string & dumper_name) { + mesh.dump(dumper_name); +} + +/* -------------------------------------------------------------------------- */ +void HeatTransferModel::dump(const std::string & dumper_name, UInt step) { + mesh.dump(dumper_name, step); +} + +/* ------------------------------------------------------------------------- */ +void HeatTransferModel::dump(const std::string & dumper_name, Real time, + UInt step) { + mesh.dump(dumper_name, time, step); +} + +/* -------------------------------------------------------------------------- */ +void HeatTransferModel::dump() { mesh.dump(); } + +/* -------------------------------------------------------------------------- */ +void HeatTransferModel::dump(UInt step) { mesh.dump(step); } + +/* -------------------------------------------------------------------------- */ +void HeatTransferModel::dump(Real time, UInt step) { mesh.dump(time, step); } + __END_AKANTU__ diff --git a/src/model/heat_transfer/heat_transfer_model.hh b/src/model/heat_transfer/heat_transfer_model.hh index c32e85634..d4e270577 100644 --- a/src/model/heat_transfer/heat_transfer_model.hh +++ b/src/model/heat_transfer/heat_transfer_model.hh @@ -1,456 +1,460 @@ /** * @file heat_transfer_model.hh * * @author Guillaume Anciaux * @author Lucas Frerot * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Tue Dec 08 2015 * * @brief Model of Heat Transfer * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_HEAT_TRANSFER_MODEL_HH__ #define __AKANTU_HEAT_TRANSFER_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" +#include "aka_memory.hh" #include "aka_types.hh" #include "aka_voigthelper.hh" -#include "aka_memory.hh" -#include "model.hh" -#include "integrator_gauss.hh" -#include "shape_lagrange.hh" #include "dumpable.hh" +#include "generalized_trapezoidal.hh" +#include "integrator_gauss.hh" +#include "model.hh" #include "parsable.hh" +#include "shape_lagrange.hh" #include "solver.hh" -#include "generalized_trapezoidal.hh" namespace akantu { - class IntegrationScheme1stOrder; +class IntegrationScheme1stOrder; } __BEGIN_AKANTU__ struct HeatTransferModelOptions : public ModelOptions { - HeatTransferModelOptions(AnalysisMethod analysis_method = _explicit_lumped_capacity ) : analysis_method(analysis_method) {} + HeatTransferModelOptions( + AnalysisMethod analysis_method = _explicit_lumped_capacity) + : analysis_method(analysis_method) {} AnalysisMethod analysis_method; }; extern const HeatTransferModelOptions default_heat_transfer_model_options; class HeatTransferModel : public Model, public DataAccessor, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - typedef FEEngineTemplate MyFEEngineType; + typedef FEEngineTemplate MyFEEngineType; - HeatTransferModel(Mesh & mesh, - UInt spatial_dimension = _all_dimensions, - const ID & id = "heat_transfer_model", - const MemoryID & memory_id = 0); + HeatTransferModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, + const ID & id = "heat_transfer_model", + const MemoryID & memory_id = 0); - virtual ~HeatTransferModel() ; + virtual ~HeatTransferModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// generic function to initialize everything ready for explicit dynamics - void initFull(const ModelOptions & options = default_heat_transfer_model_options); + void + initFull(const ModelOptions & options = default_heat_transfer_model_options); /// initialize the fem object of the boundary void initFEEngineBoundary(bool create_surface = true); /// read one material file to instantiate all the materials void readMaterials(); /// allocate all vectors void initArrays(); /// register the tags associated with the parallel synchronizer - void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL); + void initParallel(MeshPartition * partition, + DataAccessor * data_accessor = NULL); /// initialize the model void initModel(); /// init PBC synchronizer void initPBC(); /// initialize the solver and the jacobian_matrix (called by initImplicit) void initSolver(SolverOptions & solver_options); /// initialize the stuff for the implicit solver - void initImplicit(bool dynamic, SolverOptions & solver_options = _solver_no_options); + void initImplicit(bool dynamic, + SolverOptions & solver_options = _solver_no_options); /// function to print the contain of the class - virtual void printself(__attribute__ ((unused)) std::ostream & stream, - __attribute__ ((unused)) int indent = 0) const {}; + virtual void printself(__attribute__((unused)) std::ostream & stream, + __attribute__((unused)) int indent = 0) const {}; /* ------------------------------------------------------------------------ */ /* Methods for explicit */ /* ------------------------------------------------------------------------ */ public: - /// compute and get the stable time step Real getStableTimeStep(); /// compute the heat flux void updateResidual(bool compute_conductivity = false); /// calculate the lumped capacity vector for heat transfer problem void assembleCapacityLumped(); /// update the temperature from the temperature rate void explicitPred(); /// update the temperature rate from the increment void explicitCorr(); /// implicit time integration predictor void implicitPred(); /// implicit time integration corrector void implicitCorr(); - - /// solve the system in temperature rate @f$C\delta \dot T = q_{n+1} - C \dot T_{n}@f$ + /// solve the system in temperature rate @f$C\delta \dot T = q_{n+1} - C \dot + /// T_{n}@f$ /// this function needs to be run for dynamics void solveExplicitLumped(); // /// initialize the heat flux // void initializeResidual(Array &temp); // /// initialize temperature // void initializeTemperature(Array &temp); /* ------------------------------------------------------------------------ */ /* Methods for implicit */ /* ------------------------------------------------------------------------ */ public: - /** * solve Kt = q **/ void solveStatic(); /// test if the system is converged template bool testConvergence(Real tolerance, Real & error); - - /** - * solve a step (predictor + convergence loop + corrector) using the - * the given convergence method (see akantu::SolveConvergenceMethod) - * and the given convergence criteria (see - * akantu::SolveConvergenceCriteria) - **/ + /** + * solve a step (predictor + convergence loop + corrector) using the + * the given convergence method (see akantu::SolveConvergenceMethod) + * and the given convergence criteria (see + * akantu::SolveConvergenceCriteria) + **/ template bool solveStep(Real tolerance, UInt max_iteration = 100); template - bool solveStep(Real tolerance, - Real & error, - UInt max_iteration = 100, - bool do_not_factorize = false); + bool solveStep(Real tolerance, Real & error, UInt max_iteration = 100, + bool do_not_factorize = false); /// assemble the conductivity matrix void assembleConductivityMatrix(bool compute_conductivity = true); /// assemble the conductivity matrix void assembleCapacity(); /// assemble the conductivity matrix void assembleCapacity(GhostType ghost_type); /// compute the capacity on quadrature points - void computeRho(Array & rho, ElementType type, - GhostType ghost_type); - + void computeRho(Array & rho, ElementType type, GhostType ghost_type); protected: /// compute A and solve @f[ A\delta u = f_ext - f_int @f] template - void solve(Array &increment, Real block_val = 1., + void solve(Array & increment, Real block_val = 1., bool need_factorize = true, bool has_profile_changed = false); /// computation of the residual for the convergence loop void updateResidualInternal(); - -private: +private: /// compute the heat flux on ghost types - void updateResidual(const GhostType & ghost_type, bool compute_conductivity = false); + void updateResidual(const GhostType & ghost_type, + bool compute_conductivity = false); - /// calculate the lumped capacity vector for heat transfer problem (w ghosttype) + /// calculate the lumped capacity vector for heat transfer problem (w + /// ghosttype) void assembleCapacityLumped(const GhostType & ghost_type); /// assemble the conductivity matrix (w/ ghost type) template void assembleConductivityMatrix(const GhostType & ghost_type, - bool compute_conductivity = true); + bool compute_conductivity = true); /// assemble the conductivity matrix template void assembleConductivityMatrix(const ElementType & type, - const GhostType & ghost_type, - bool compute_conductivity = true); + const GhostType & ghost_type, + bool compute_conductivity = true); /// compute the conductivity tensor for each quadrature point in an array void computeConductivityOnQuadPoints(const GhostType & ghost_type); /// compute vector k \grad T for each quadrature point - void computeKgradT(const GhostType & ghost_type,bool compute_conductivity); + void computeKgradT(const GhostType & ghost_type, bool compute_conductivity); /// compute the thermal energy Real computeThermalEnergyByNode(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbDataForElements(const Array & elements, - SynchronizationTag tag) const; + SynchronizationTag tag) const; inline void packElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag) const; + const Array & elements, + SynchronizationTag tag) const; inline void unpackElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag); + const Array & elements, + SynchronizationTag tag); inline UInt getNbDataToPack(SynchronizationTag tag) const; inline UInt getNbDataToUnpack(SynchronizationTag tag) const; - inline void packData(CommunicationBuffer & buffer, - const UInt index, - SynchronizationTag tag) const; - inline void unpackData(CommunicationBuffer & buffer, - const UInt index, - SynchronizationTag tag); + inline void packData(CommunicationBuffer & buffer, const UInt index, + SynchronizationTag tag) const; + inline void unpackData(CommunicationBuffer & buffer, const UInt index, + SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: - virtual dumper::Field * createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag); + const std::string & group_name, + bool padding_flag); virtual dumper::Field * createNodalFieldBool(const std::string & field_name, - const std::string & group_name, - bool padding_flag); + const std::string & group_name, + bool padding_flag); + + virtual dumper::Field * createElementalField(const std::string & field_name, + const std::string & group_name, + bool padding_flag, + const UInt & spatial_dimension, + const ElementKind & kind); + + virtual void dump(const std::string & dumper_name); + + virtual void dump(const std::string & dumper_name, UInt step); + + virtual void dump(const std::string & dumper_name, Real time, UInt step); + + virtual void dump(); + virtual void dump(UInt step); - virtual dumper::Field * createElementalField(const std::string & field_name, - const std::string & group_name, - bool padding_flag, - const UInt & spatial_dimension, - const ElementKind & kind); + virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - inline FEEngine & getFEEngineBoundary(const std::string & name = ""); AKANTU_GET_MACRO(Density, density, Real); AKANTU_GET_MACRO(Capacity, capacity, Real); /// get the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step AKANTU_SET_MACRO(TimeStep, time_step, Real); /// get the assembled heat flux - AKANTU_GET_MACRO(Residual, *residual, Array&); + AKANTU_GET_MACRO(Residual, *residual, Array &); /// get the lumped capacity - AKANTU_GET_MACRO(CapacityLumped, *capacity_lumped, Array&); + AKANTU_GET_MACRO(CapacityLumped, *capacity_lumped, Array &); /// get the boundary vector - AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array&); + AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get conductivity matrix - AKANTU_GET_MACRO(ConductivityMatrix, *conductivity_matrix, const SparseMatrix&); + AKANTU_GET_MACRO(ConductivityMatrix, *conductivity_matrix, + const SparseMatrix &); /// get the external heat rate vector - AKANTU_GET_MACRO(ExternalHeatRate, *external_heat_rate, Array&); + AKANTU_GET_MACRO(ExternalHeatRate, *external_heat_rate, Array &); /// get the temperature gradient - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureGradient, temperature_gradient, Real); + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureGradient, + temperature_gradient, Real); /// get the conductivity on q points - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConductivityOnQpoints, conductivity_on_qpoints, Real); + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConductivityOnQpoints, + conductivity_on_qpoints, Real); /// get the conductivity on q points - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureOnQpoints, temperature_on_qpoints, Real); + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureOnQpoints, + temperature_on_qpoints, Real); /// internal variables - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(KGradtOnQpoints, k_gradt_on_qpoints, Real); + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(KGradtOnQpoints, k_gradt_on_qpoints, + Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(IntBtKgT, int_bt_k_gT, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(BtKgT, bt_k_gT, Real); /// get the temperature AKANTU_GET_MACRO(Temperature, *temperature, Array &); /// get the temperature derivative AKANTU_GET_MACRO(TemperatureRate, *temperature_rate, Array &); /// get the equation number Array AKANTU_GET_MACRO(EquationNumber, *equation_number, const Array &); /// get the energy denominated by thermal - Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); + Real getEnergy(const std::string & energy_id, const ElementType & type, + UInt index); /// get the energy denominated by thermal Real getEnergy(const std::string & energy_id); /// get the thermal energy for a given element Real getThermalEnergy(const ElementType & type, UInt index); /// get the thermal energy for a given element Real getThermalEnergy(); protected: /* ----------------------------------------------------------------------- */ - template - void getThermalEnergy(iterator Eth, - Array::const_iterator T_it, - Array::const_iterator T_end) const; + template + void getThermalEnergy(iterator Eth, Array::const_iterator T_it, + Array::const_iterator T_end) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: - - /// number of iterations + /// number of iterations UInt n_iter; IntegrationScheme1stOrder * integrator; /// time step Real time_step; /// temperatures array Array * temperature; /// temperatures derivatives array Array * temperature_rate; /// increment array (@f$\delta \dot T@f$ or @f$\delta T@f$) Array * increment; /// conductivity matrix SparseMatrix * conductivity_matrix; /// capacity matrix - SparseMatrix *capacity_matrix; + SparseMatrix * capacity_matrix; /// jacobian matrix SparseMatrix * jacobian_matrix; /// the density Real density; /// the speed of the changing temperature ElementTypeMapArray temperature_gradient; /// temperature field on quadrature points ElementTypeMapArray temperature_on_qpoints; /// conductivity tensor on quadrature points ElementTypeMapArray conductivity_on_qpoints; /// vector k \grad T on quad points ElementTypeMapArray k_gradt_on_qpoints; /// vector \int \grad N k \grad T ElementTypeMapArray int_bt_k_gT; /// vector \grad N k \grad T ElementTypeMapArray bt_k_gT; /// external flux vector Array * external_heat_rate; /// residuals array Array * residual; /// position of a dof in the K matrix Array * equation_number; - //lumped vector + // lumped vector Array * capacity_lumped; /// boundary vector Array * blocked_dofs; - //realtime + // realtime Real time; - ///capacity + /// capacity Real capacity; - //conductivity matrix + // conductivity matrix Matrix conductivity; - //linear variation of the conductivity (for temperature dependent conductivity) + // linear variation of the conductivity (for temperature dependent + // conductivity) Real conductivity_variation; // reference temperature for the interpretation of temperature variation Real T_ref; - //the biggest parameter of conductivity matrix + // the biggest parameter of conductivity matrix Real conductivitymax; /// thermal energy by element ElementTypeMapArray thermal_energy; /// Solver Solver * solver; /// analysis method AnalysisMethod method; /// pointer to the pbc synchronizer PBCSynchronizer * pbc_synch; - }; - /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ -#if defined (AKANTU_INCLUDE_INLINE_IMPL) -# include "heat_transfer_model_inline_impl.cc" +#if defined(AKANTU_INCLUDE_INLINE_IMPL) +#include "heat_transfer_model_inline_impl.cc" #endif /// standard output stream operator -inline std::ostream & operator <<(std::ostream & stream, const HeatTransferModel & _this) -{ +inline std::ostream & operator<<(std::ostream & stream, + const HeatTransferModel & _this) { _this.printself(stream); return stream; } - __END_AKANTU__ - - #endif /* __AKANTU_HEAT_TRANSFER_MODEL_HH__ */