diff --git a/packages/heat_transfer.cmake b/packages/heat_transfer.cmake index c414654c6..3b6e2052f 100644 --- a/packages/heat_transfer.cmake +++ b/packages/heat_transfer.cmake @@ -1,43 +1,44 @@ #=============================================================================== # @file heat_transfer.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Mon Mar 30 2015 # # @brief package description for heat transfer # # # @section LICENSE # # Copyright (©) 2010-2021 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 . # #=============================================================================== package_declare(heat_transfer DESCRIPTION "Use Heat Transfer package of Akantu") package_declare_sources(heat_transfer model/heat_transfer/heat_transfer_model.cc model/heat_transfer/heat_transfer_model.hh model/heat_transfer/heat_transfer_model_inline_impl.hh model/heat_transfer_interface/heat_transfer_interface_model.cc model/heat_transfer_interface/heat_transfer_interface_model.hh - ) + model/heat_transfer_interface/heat_transfer_interface_model_inline_impl.hh +) diff --git a/src/model/heat_transfer_interface/heat_transfer_interface_model.cc b/src/model/heat_transfer_interface/heat_transfer_interface_model.cc index 16660a971..a7b42afa2 100644 --- a/src/model/heat_transfer_interface/heat_transfer_interface_model.cc +++ b/src/model/heat_transfer_interface/heat_transfer_interface_model.cc @@ -1,805 +1,822 @@ /** * @file heat_transfer_interface_model.hh * * @author Emil Gallyamov * * @date creation: Thu Apr 13 2023 * @date last modification: Thu Apr 13 2023 * * @brief Implementation of HeatTransferInterfaceModel class * * * @section LICENSE * * Copyright (©) 2010-2021 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 "heat_transfer_interface_model.hh" #include "dumpable_inline_impl.hh" #include "element_class.hh" #include "element_synchronizer.hh" #include "fe_engine_template.hh" #include "generalized_trapezoidal.hh" #include "group_manager_inline_impl.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "parser.hh" #include "shape_cohesive.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ namespace akantu { namespace heat_transfer_interface { namespace details { class ComputeRhoFunctor { public: ComputeRhoFunctor(const HeatTransferInterfaceModel & model) : model(model){}; void operator()(Matrix & rho, const Element & el) { auto opening_it = model.getOpening(el.type, el.ghost_type).begin(); rho.set(model.getCapacityInCrack() * model.getDensityInCrack() * opening_it[el.element]); } private: const HeatTransferInterfaceModel & model; }; } // namespace details } // namespace heat_transfer_interface /* -------------------------------------------------------------------------- */ HeatTransferInterfaceModel::HeatTransferInterfaceModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager) : HeatTransferModel(mesh, dim, id, ModelType::_heat_transfer_interface_model, dof_manager), opening_on_qpoints("opening_on_qpoints", id), opening_rate("opening_rate", id), k_long_w("k_long_w", id), k_perp_over_w("k_perp_over_w", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("InterfacesFEEngine", mesh, Model::spatial_dimension); this->mesh.registerDumper("heat_interfaces", id); this->mesh.addDumpMeshToDumper("heat_interfaces", mesh, spatial_dimension, _not_ghost, _ek_cohesive); if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_htm_gradient_temperature); + + auto & mesh_facets = this->mesh.getMeshFacets(); + auto & facet_synchronizer = mesh_facets.getElementSynchronizer(); + const auto & cfacet_synchronizer = facet_synchronizer; + // update the cohesive element synchronizer + cohesive_synchronizer->updateSchemes([&](auto && scheme, auto && proc, + auto && direction) { + auto & facet_scheme = + cfacet_synchronizer.getCommunications().getScheme(proc, direction); + + for (auto && facet : facet_scheme) { + const auto & cohesive_element = const_cast(mesh_facets) + .getElementToSubelement(facet)[1]; + + if (cohesive_element == ElementNull or + cohesive_element.kind() != _ek_cohesive) { + continue; + } + scheme.push_back(cohesive_element); + } + }); } + this->registerParam("longitudinal_conductivity", longitudinal_conductivity, _pat_parsmod); this->registerParam("transversal_conductivity", transversal_conductivity, _pat_parsmod); this->registerParam("default_opening", default_opening, 1.e-5, _pat_parsmod); this->registerParam("capacity_in_crack", capacity_in_crack, _pat_parsmod); this->registerParam("density_in_crack", density_in_crack, _pat_parsmod); this->registerParam("use_opening_rate", use_opening_rate, _pat_parsmod); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::initModel() { Parent::initModel(); auto & fem = this->getFEEngine("InterfacesFEEngine"); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); this->temperature_gradient.initialize(fem, _nb_component = spatial_dimension); this->k_long_w.initialize(fem, _nb_component = (spatial_dimension - 1) * (spatial_dimension - 1)); this->k_perp_over_w.initialize(fem, _nb_component = 1); this->opening_on_qpoints.initialize(fem, _nb_component = 1); if (use_opening_rate) { opening_rate.initialize(fem, _nb_component = 1); } } /* -------------------------------------------------------------------------- */ HeatTransferInterfaceModel::~HeatTransferInterfaceModel() = default; /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::predictor() { Parent::predictor(); /// to force HeatTransferModel recompute conductivity matrix ++this->conductivity_release[_not_ghost]; /// same for HeatTransferInterfaceModel ++opening_release; } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::readMaterials() { auto sect = this->getParserSection(); if (not std::get<1>(sect)) { const auto & section = std::get<0>(sect); this->parseSection(section); } this->conductivity_on_qpoints.set(conductivity); this->initial_conductivity_array.set(conductivity); this->density_array.set(density); this->capacity_array.set(capacity); this->conductivity_variation_array.set(conductivity_variation); opening_on_qpoints.set(default_opening); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::assembleCapacity() { AKANTU_DEBUG_IN(); auto ghost_type = _not_ghost; Parent::assembleCapacity(); auto & fem_interface = this->getFEEngine("InterfacesFEEngine"); heat_transfer_interface::details::ComputeRhoFunctor rho_functor(*this); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { fem_interface.assembleFieldMatrix(rho_functor, "M", "temperature", this->getDOFManager(), type, ghost_type); } need_to_reassemble_capacity = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void HeatTransferInterfaceModel::assembleCapacityLumped(GhostType ghost_type) // { // AKANTU_DEBUG_IN(); // Parent::assembleCapacityLumped(ghost_type); // auto & fem_interface = // this->getFEEngineClass("InterfacesFEEngine"); // heat_transfer_interface::details::ComputeRhoFunctor compute_rho(*this); // for (auto && type : // mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { // fem_interface.assembleFieldLumped(compute_rho, "M", "temperature", // this->getDOFManager(), type, // ghost_type); // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::assembleConductivityMatrix() { AKANTU_DEBUG_IN(); Parent::assembleConductivityMatrix(); this->computeKLongOnQuadPoints(_not_ghost); this->computeKTransOnQuadPoints(_not_ghost); if ((long_conductivity_release[_not_ghost] == this->crack_conductivity_matrix_release) and (perp_conductivity_release[_not_ghost] == this->crack_conductivity_matrix_release)) { return; } auto & fem_interface = this->getFEEngine("InterfacesFEEngine"); for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_cohesive)) { auto nb_element = mesh.getNbElement(type); if (nb_element == 0) continue; auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem_interface.getNbIntegrationPoints(type); auto && shape_derivatives = fem_interface.getShapesDerivatives(type); auto && shapes = fem_interface.getShapes(type); auto A = ExtendingOperators::getAveragingOperator(type); auto C = ExtendingOperators::getDifferencingOperator(type); // tangent_conductivity_matrix.zero(); auto && long_cond = this->k_long_w(type); auto && perp_cond = this->k_perp_over_w(type); Array at_bt_k_long_b_a(nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "A^t*B^t*k_long*B*A"); Array ct_nt_k_perp_n_c(nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "C^t*N^t*k_perp*N*C"); Array tangent_conductivity( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "tangent_conductivity"); Matrix B_A(spatial_dimension - 1, nb_nodes_per_element); Matrix k_long_B_A(spatial_dimension - 1, nb_nodes_per_element); // Matrix N(spatial_dimension, nb_nodes_per_element); Matrix N_C(1, nb_nodes_per_element); for (auto && data : zip(make_view(long_cond, spatial_dimension - 1, spatial_dimension - 1), make_view(shape_derivatives, spatial_dimension - 1, nb_nodes_per_element / 2), make_view(at_bt_k_long_b_a, nb_nodes_per_element, nb_nodes_per_element), perp_cond, make_view(shapes, 1, nb_nodes_per_element / 2), make_view(ct_nt_k_perp_n_c, nb_nodes_per_element, nb_nodes_per_element), make_view(tangent_conductivity, nb_nodes_per_element, nb_nodes_per_element))) { // assemble conductivity contribution along crack const auto & k_long = std::get<0>(data); const auto & B = std::get<1>(data); auto & At_Bt_k_long_B_A = std::get<2>(data); B_A.mul(B, A); k_long_B_A.mul(k_long, B_A); At_Bt_k_long_B_A.mul(k_long_B_A, B_A); // assemble conductivity contribution perpendicular to the crack const auto & k_perp = std::get<3>(data); const auto & N = std::get<4>(data); auto & Ct_Nt_k_perp_N_C = std::get<5>(data); auto & tangent = std::get<6>(data); N_C.mul(N, C); Ct_Nt_k_perp_N_C.mul(N_C, N_C, k_perp); // summing two contributions of tangent operator tangent = At_Bt_k_long_B_A + Ct_Nt_k_perp_N_C; } auto K_e = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_e"); fem_interface.integrate(tangent_conductivity, *K_e, nb_nodes_per_element * nb_nodes_per_element, type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "temperature", *K_e, type, _not_ghost, _unsymmetric); } this->crack_conductivity_matrix_release = long_conductivity_release[_not_ghost]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::computeGradAndDeltaT(GhostType ghost_type) { auto & fem_interface = this->getFEEngineClass("InterfacesFEEngine"); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto nb_quadrature_points = fem_interface.getNbIntegrationPoints(type); auto nb_element = mesh.getNbElement(type, ghost_type); auto & gradient = temperature_gradient(type, ghost_type); auto surface_gradient = std::make_unique>( nb_element * nb_quadrature_points, spatial_dimension - 1, "surface_gradient"); auto jump = std::make_unique>(nb_element * nb_quadrature_points, 1, "temperature_jump"); this->getFEEngine("InterfacesFEEngine") .gradientOnIntegrationPoints(*temperature, *surface_gradient, 1, type, ghost_type); #define COMPUTE_JUMP(type) \ fem_interface.getShapeFunctions() \ .interpolateOnIntegrationPoints( \ *temperature, *jump, 1, ghost_type); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_JUMP); #undef COMPUTE_JUMP for (auto && data : zip(make_view(gradient, spatial_dimension), make_view(*surface_gradient, spatial_dimension - 1), *jump)) { auto & grad = std::get<0>(data); auto & surf_grad = std::get<1>(data); auto & delta = std::get<2>(data); for (auto i : arange(spatial_dimension - 1)) { grad(i) = surf_grad(i); } grad(spatial_dimension - 1) = delta; } } } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::computeKLongOnQuadPoints( GhostType ghost_type) { // if opening did not change, longitudinal conductivity will neither if (opening_release == long_conductivity_release[ghost_type]) { return; } for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto & opening = opening_on_qpoints(type, ghost_type); auto & long_cond_w = k_long_w(type, ghost_type); UInt dim = spatial_dimension - 1; Matrix long_cond_vect(dim, dim); for (auto i : arange(dim)) { long_cond_vect(i, i) = this->longitudinal_conductivity; } for (auto && tuple : zip(make_view(long_cond_w, dim, dim), opening)) { auto & k = std::get<0>(tuple); auto & w = std::get<1>(tuple); k = long_cond_vect * w; } } long_conductivity_release[ghost_type] = opening_release; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::computeKTransOnQuadPoints( GhostType ghost_type) { // if opening did not change, longitudinal conductivity will neither if (opening_release == perp_conductivity_release[ghost_type]) { return; } for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto & opening = opening_on_qpoints(type, ghost_type); auto & trans_cond_over_w = k_perp_over_w(type, ghost_type); for (auto && tuple : zip(trans_cond_over_w, opening)) { auto & k_perp = std::get<0>(tuple); auto & w = std::get<1>(tuple); k_perp = transversal_conductivity / w; } } perp_conductivity_release[ghost_type] = opening_release; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::assembleInternalHeatRate() { AKANTU_DEBUG_IN(); Parent::assembleInternalHeatRate(); computeGradAndDeltaT(_not_ghost); // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_htm_gradient_temperature); computeLongHeatRate(_not_ghost); computeTransHeatRate(_not_ghost); computeInertialHeatRate(_not_ghost); // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(SynchronizationTag::_htm_gradient_temperature); computeLongHeatRate(_ghost); computeTransHeatRate(_ghost); computeInertialHeatRate(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::computeLongHeatRate(GhostType ghost_type) { AKANTU_DEBUG_IN(); computeKLongOnQuadPoints(ghost_type); auto & internal_rate = const_cast &>(this->getInternalHeatRate()); auto & fem = this->getFEEngine("InterfacesFEEngine"); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { const auto & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); auto & gradient = temperature_gradient(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto A = ExtendingOperators::getAveragingOperator(type); Matrix B_A(spatial_dimension - 1, nb_nodes_per_element); auto k_long_w_gradT = std::make_unique>(nb_element * nb_quadrature_points, spatial_dimension - 1, "k_long_grad_T"); auto bt_k_w_gradT = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element); for (auto && values : zip(make_view(k_long_w(type, ghost_type), spatial_dimension - 1, spatial_dimension - 1), make_view(gradient, spatial_dimension), make_view(*k_long_w_gradT, spatial_dimension - 1), make_view(shapes_derivatives, spatial_dimension - 1, nb_nodes_per_element / 2), make_view(*bt_k_w_gradT, nb_nodes_per_element))) { const auto & k_w = std::get<0>(values); const auto & full_gradT = std::get<1>(values); Vector surf_gradT(spatial_dimension - 1); - bool zero_grad = true; for (auto i : arange(spatial_dimension - 1)) { surf_gradT(i) = full_gradT(i); - if (full_gradT(i) != 0) - zero_grad = false; - } - if (zero_grad and ghost_type == _ghost) { - std::cout << "zero grad" << std::endl; } + auto & k_w_gradT = std::get<2>(values); k_w_gradT.mul(k_w, surf_gradT); /// compute @f$B_a t_i@f$ const auto & B = std::get<3>(values); auto & At_Bt_vector = std::get<4>(values); B_A.mul(B, A); At_Bt_vector.template mul(B_A, k_w_gradT); } auto long_heat_rate = std::make_unique>( nb_element, nb_nodes_per_element, "long_heat_rate"); fem.integrate(*bt_k_w_gradT, *long_heat_rate, nb_nodes_per_element, type, ghost_type); /// assemble this->getDOFManager().assembleElementalArrayLocalArray( *long_heat_rate, internal_rate, type, ghost_type, 1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::computeTransHeatRate(GhostType ghost_type) { AKANTU_DEBUG_IN(); computeKTransOnQuadPoints(ghost_type); auto & internal_rate = const_cast &>(this->getInternalHeatRate()); auto & fem = this->getFEEngine("InterfacesFEEngine"); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { const auto & shapes = fem.getShapes(type, ghost_type); auto size_of_shapes = shapes.getNbComponent(); auto & gradient = temperature_gradient(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto natural_dimension = spatial_dimension - 1; // difference computing operator auto C = ExtendingOperators::getDifferencingOperator(type); auto nt_k_deltaT = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element); auto k_perp_deltaT_over_w = std::make_unique>( nb_element * nb_quadrature_points, 1, "k_perp_deltaT"); for (auto && values : zip(*k_perp_deltaT_over_w, k_perp_over_w(type, ghost_type), make_view(gradient, spatial_dimension), make_view(shapes, size_of_shapes), make_view(*nt_k_deltaT, nb_nodes_per_element))) { auto & _k_perp_deltaT_over_w = std::get<0>(values); auto & _k_perp_over_w = std::get<1>(values); const auto & full_gradT = std::get<2>(values); const auto & N = std::get<3>(values); auto & Nt_k_deltaT = std::get<4>(values); _k_perp_deltaT_over_w = _k_perp_over_w * full_gradT(spatial_dimension - 1); Nt_k_deltaT.mul(C, N, _k_perp_deltaT_over_w); } auto perp_heat_rate = std::make_unique>( nb_element, nb_nodes_per_element, "perp_heat_rate"); fem.integrate(*nt_k_deltaT, *perp_heat_rate, nb_nodes_per_element, type, ghost_type); /// assemble this->getDOFManager().assembleElementalArrayLocalArray( *perp_heat_rate, internal_rate, type, ghost_type, 1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------*/ void HeatTransferInterfaceModel::computeInertialHeatRate(GhostType ghost_type) { AKANTU_DEBUG_IN(); if (not this->use_opening_rate) return; auto & internal_rate = const_cast &>(this->getInternalHeatRate()); auto & fem_interface = this->getFEEngineClass("InterfacesFEEngine"); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { const auto & shapes = fem_interface.getShapes(type, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem_interface.getNbIntegrationPoints(type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto natural_dimension = spatial_dimension - 1; // averaging operator auto A = ExtendingOperators::getAveragingOperator(type); auto && opening_rates = this->opening_rate(type, ghost_type); auto * nt_opening_rate = new Array(nb_element * nb_quadrature_points, nb_nodes_per_element); for (auto && values : zip(make_view(shapes, size_of_shapes), make_view(*nt_opening_rate, nb_nodes_per_element), opening_rates)) { const auto & N = std::get<0>(values); auto & Nt_opening_rate = std::get<1>(values); auto & open_rate = std::get<2>(values); Nt_opening_rate.mul(A, N, open_rate); } auto * inertial_heat_rate = new Array(nb_element, nb_nodes_per_element, "inertial_heat_rate"); fem_interface.integrate(*nt_opening_rate, *inertial_heat_rate, nb_nodes_per_element, type, ghost_type); delete nt_opening_rate; /// assemble this->getDOFManager().assembleElementalArrayLocalArray( *inertial_heat_rate, internal_rate, type, ghost_type, 1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const ShapeLagrange<_ek_cohesive> & HeatTransferInterfaceModel::getShapeFunctionsCohesive() { return this->getFEEngineClass("InterfacesFEEngine") .getShapeFunctions(); } /* -------------------------------------------------------------------------- */ Real HeatTransferInterfaceModel::getStableTimeStep() { AKANTU_DEBUG_IN(); auto HTM_time_step = Parent::getStableTimeStep(); Real el_size; Real min_el_size = std::numeric_limits::max(); Real max_opening_on_qpoint = std::numeric_limits::min(); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_cohesive)) { auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); auto & opening = opening_on_qpoints(type, _not_ghost); auto nb_quad_per_element = opening.getNbComponent(); auto mesh_dim = this->mesh.getSpatialDimension(); const auto facet_type = Mesh::getFacetType(type); Array coord(0, nb_nodes_per_element * mesh_dim); FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, _not_ghost); for (auto && data : zip(make_view(coord, mesh_dim, nb_nodes_per_element), make_view(opening, nb_quad_per_element))) { Matrix & el_coord = std::get<0>(data); Vector & el_opening = std::get<1>(data); el_size = getFEEngine().getElementInradius(el_coord, facet_type); min_el_size = std::min(min_el_size, el_size); max_opening_on_qpoint = std::max(max_opening_on_qpoint, el_opening.norm()); } AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size << " and the max opening is : " << max_opening_on_qpoint); } Real min_dt = 2. * min_el_size * min_el_size / 4 * density_in_crack * capacity_in_crack / longitudinal_conductivity; mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); min_dt = std::min(min_dt, HTM_time_step); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::setTimeStep(Real time_step, const ID & solver_id) { Parent::setTimeStep(time_step, solver_id); this->mesh.getDumper("heat_interfaces").setTimeStep(time_step); } // /* // -------------------------------------------------------------------------- // */ void HeatTransferInterfaceModel::assignPropertyToPhysicalGroup( // const std::string & property_name, const std::string & group_name, // Real value) { // AKANTU_DEBUG_ASSERT(property_name != "conductivity", // "Scalar was provided instead of a conductivity // matrix"); // auto && el_group = mesh.getElementGroup(group_name); // auto & fem = this->getFEEngine(); // for (auto ghost_type : ghost_types) { // for (auto && type : // mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { // auto nb_quadrature_points = fem.getNbIntegrationPoints(type); // auto && elements = el_group.getElements(type, ghost_type); // Array::scalar_iterator field_it; // if (property_name == "density") { // field_it = density_array(type, ghost_type).begin(); // } else if (property_name == "capacity") { // field_it = capacity_array(type, ghost_type).begin(); // } else { // AKANTU_EXCEPTION(property_name + // " is not a valid material property name."); // } // for (auto && el : elements) { // for (auto && qpoint : arange(nb_quadrature_points)) { // field_it[el * nb_quadrature_points + qpoint] = value; // } // } // need_to_reassemble_capacity = true; // need_to_reassemble_capacity_lumped = true; // } // } // } // /* // -------------------------------------------------------------------------- // */ void HeatTransferInterfaceModel::assignPropertyToPhysicalGroup( // const std::string & property_name, const std::string & group_name, // Matrix cond_matrix) { // AKANTU_DEBUG_ASSERT(property_name == "conductivity", // "When updating material parameters, only conductivity // " "accepts matrix as an input"); // auto && el_group = mesh.getElementGroup(group_name); // auto & fem = this->getFEEngine(); // auto dim = this->getMesh().getSpatialDimension(); // for (auto ghost_type : ghost_types) { // for (auto && type : // mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { // auto nb_quadrature_points = fem.getNbIntegrationPoints(type); // auto && elements = el_group.getElements(type, ghost_type); // auto init_cond_it = // make_view(initial_conductivity_array(type, ghost_type), dim, dim) // .begin(); // auto cond_on_quad_it = // make_view(conductivity_on_qpoints(type, ghost_type), dim, dim) // .begin(); // for (auto && el : elements) { // for (auto && qpoint : arange(nb_quadrature_points)) { // init_cond_it[el * nb_quadrature_points + qpoint] = cond_matrix; // cond_on_quad_it[el * nb_quadrature_points + qpoint] = // cond_matrix; // } // } // } // conductivity_release[ghost_type] += 1; // } // } /* -------------------------------------------------------------------------*/ std::shared_ptr HeatTransferInterfaceModel::createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind element_kind) { std::shared_ptr field; auto getNbDataPerElem = [&](auto & field) { const auto & fe_engine = getFEEngine("InterfacesFEEngine"); ElementTypeMap res; for (auto ghost_type : ghost_types) { for (auto & type : field.elementTypes(ghost_type)) { auto nb_quadrature_points = fe_engine.getNbIntegrationPoints(type, ghost_type); res(type, ghost_type) = field(type, ghost_type).getNbComponent() * nb_quadrature_points; } } return res; }; if (element_kind == _ek_regular) { field = Parent::createElementalField(field_name, group_name, padding_flag, spatial_dimension, element_kind); } else if (element_kind == _ek_cohesive) { if (field_name == "partitions") { field = mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); } else if (field_name == "opening_on_qpoints") { auto nb_data_per_elem = getNbDataPerElem(opening_on_qpoints); field = mesh.createElementalField( opening_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "temperature_gradient") { auto nb_data_per_elem = getNbDataPerElem(temperature_gradient); field = mesh.createElementalField( temperature_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } } return field; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/heat_transfer_interface/heat_transfer_interface_model.hh b/src/model/heat_transfer_interface/heat_transfer_interface_model.hh index b4d14affc..a903d54d0 100644 --- a/src/model/heat_transfer_interface/heat_transfer_interface_model.hh +++ b/src/model/heat_transfer_interface/heat_transfer_interface_model.hh @@ -1,253 +1,245 @@ /** * @file heat_transfer_interface_model.hh * * @author Emil Gallyamov * * @date creation: Thu Apr 13 2023 * @date last modification: Thu Apr 13 2023 * * @brief Model of Heat Transfer expanded to heat diffusion along interfaces * * * @section LICENSE * * Copyright (©) 2010-2021 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 "heat_transfer_model.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_HEAT_TRANSFER_INTERFACE_MODEL_HH_ #define AKANTU_HEAT_TRANSFER_INTERFACE_MODEL_HH_ namespace akantu { class HeatTransferInterfaceModel : public HeatTransferModel { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineCohesiveType = FEEngineTemplate; using Parent = HeatTransferModel; HeatTransferInterfaceModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "heat_transfer_interface_model", std::shared_ptr dof_manager = nullptr); ~HeatTransferInterfaceModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// read one material file to instantiate all the materials void readMaterials() override; /// initialize the model void initModel() override; void predictor() override; /* ------------------------------------------------------------------------ */ /* Methods for explicit */ /* ------------------------------------------------------------------------ */ public: /// compute and get the stable time step Real getStableTimeStep() override; /// set the stable timestep void setTimeStep(Real time_step, const ID & solver_id = "") override; public: // /// calculate the lumped capacity vector for heat transfer problem // void assembleCapacityLumped() override; public: /// compute the internal heat flux void assembleInternalHeatRate() override; /// assemble the conductivity matrix void assembleConductivityMatrix() override; // /// compute normals on cohesive elements // void computeNormal(const Array & position, Array & normal, // ElementType type, GhostType ghost_type = _not_ghost); // /// compute basis alliged with cohesive elements // void computeBasis(const Array & position, Array & basis, // ElementType type, GhostType ghost_type = _not_ghost); /// assemble the capacity matrix void assembleCapacity() override; /// compute the capacity on quadrature points void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// get FEEngine for cohesive elements const ShapeLagrange<_ek_cohesive> & getShapeFunctionsCohesive(); public: // /// asign material properties to physical groups // void assignPropertyToPhysicalGroup(const std::string & property_name, // const std::string & group_name, // Real value); // /// asign material properties to physical groups // void assignPropertyToPhysicalGroup(const std::string & property_name, // const std::string & group_name, // Matrix cond_matrix); private: /// compute the integrated longitudinal conductivity matrix (or scalar in /// 2D) for each quadrature point void computeKLongOnQuadPoints(GhostType ghost_type); /// compute the transversal conductivity scalar devided by opening void computeKTransOnQuadPoints(GhostType ghost_type); void computeGradAndDeltaT(GhostType ghost_type); /// compute internal heat rate along the crack and assemble it into internal /// forces void computeLongHeatRate(GhostType ghost_type); /// compute internal heat rate perpendicular to the crack and assemble it into /// internal forces void computeTransHeatRate(GhostType ghost_type); /// compute heat rate contributions due to the opening/closure rate /// (dw/dt) void computeInertialHeatRate(GhostType ghost_type); /// compute an averaging operator of size (nb_nodes_per_itp_type x /// nb_nodes_per_cohesive_type) Matrix getAveragingOperator(const ElementType & type, const UInt & nb_degree_of_freedom = 1); /// compute an differencing operator of size (nb_nodes_per_itp_type x /// nb_nodes_per_cohesive_type) Matrix getDifferencingOperator(const ElementType & type, const UInt & nb_degree_of_freedom = 1); // /// compute the thermal energy // Real computeThermalEnergyByNode(); protected: /// calculate the lumped capacity vector for heat transfer problem (w /// ghost type) void assembleCapacityLumped(GhostType ghost_type) override { AKANTU_TO_IMPLEMENT(); }; /* ------------------------------------------------------------------------ */ + /* Data Accessor inherited members */ + /* ------------------------------------------------------------------------ */ +public: + inline UInt getNbData(const Array & elements, + const SynchronizationTag & tag) const override; + inline void packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const override; + inline void unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) override; + /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(TransversalConductivity, transversal_conductivity, Real); AKANTU_GET_MACRO(LongitudinalConductivity, longitudinal_conductivity, Real); AKANTU_GET_MACRO(CapacityInCrack, capacity_in_crack, Real); AKANTU_GET_MACRO(DensityInCrack, density_in_crack, Real); AKANTU_GET_MACRO(DefaultOpening, default_opening, Real); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TransversalConductivityOnQpoints, k_perp_over_w, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(LongitudinalConductivityOnQpoints, k_long_w, Real); /// get the aperture on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening_on_qpoints, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Opening, opening_on_qpoints, Real); -public: - inline UInt getNbData(const Array & elements, - const SynchronizationTag & tag) const override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: Real transversal_conductivity; Real longitudinal_conductivity; Real default_opening; Real capacity_in_crack; Real density_in_crack; /// crack opening interpolated on quadrature points ElementTypeMapArray opening_on_qpoints; /// opening rate at quadrature points ElementTypeMapArray opening_rate; /// longitudinal conductivity tensor (or a scalar in 2D) multiplied by opening /// on quadrature points ElementTypeMapArray k_long_w; /// transversal conductivity scalar devided by opening on quadrature points ElementTypeMapArray k_perp_over_w; /// @brief boolean enabling opening-rate term into internal heat rate bool use_opening_rate{false}; std::unordered_map long_conductivity_release{{_not_ghost, 0}, {_ghost, 0}}; std::unordered_map perp_conductivity_release{{_not_ghost, 0}, {_ghost, 0}}; UInt crack_conductivity_matrix_release{UInt(-1)}; UInt opening_release{0}; /// cohesive elements synchronizer std::unique_ptr cohesive_synchronizer; }; +} // namespace akantu /* -------------------------------------------------------------------------- */ -inline UInt -HeatTransferInterfaceModel::getNbData(const Array & elements, - const SynchronizationTag & tag) const { - AKANTU_DEBUG_IN(); - - UInt size = 0; - UInt nb_nodes_per_element = 0; - Array::const_iterator it = elements.begin(); - Array::const_iterator end = elements.end(); - for (; it != end; ++it) { - const Element & el = *it; - if (el.kind() == _ek_cohesive) - std::cout << "Cohesive syncronized" << std::endl; - nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); - } - - size += Parent::getNbData(elements, tag); - AKANTU_DEBUG_OUT(); - return size; -} -} // namespace akantu +/* inline functions */ +/* -------------------------------------------------------------------------- */ +#include "heat_transfer_interface_model_inline_impl.hh" #endif /* AKANTU_HEAT_TRANSFER_INTERFACE_MODEL_HH_ */ diff --git a/src/model/heat_transfer_interface/heat_transfer_interface_model_inline_impl.hh b/src/model/heat_transfer_interface/heat_transfer_interface_model_inline_impl.hh new file mode 100644 index 000000000..ddc391769 --- /dev/null +++ b/src/model/heat_transfer_interface/heat_transfer_interface_model_inline_impl.hh @@ -0,0 +1,141 @@ +/** + * @file heat_transfer_interface_model_inline_impl.hh + * + * @author Emil Gallyamov + * + * @date creation: Wed May 24 2023 + * @date last modification: Wed May 24 2023 + * + * @brief Implementation of the inline functions of the + * HeatTransferInterfaceModel class + * + * + * @section LICENSE + * + * Copyright (©) 2010-2021 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_INTERFACE_MODEL_INLINE_IMPL_HH_ +#define AKANTU_HEAT_TRANSFER_INTERFACE_MODEL_INLINE_IMPL_HH_ + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +inline UInt +HeatTransferInterfaceModel::getNbData(const Array & elements, + const SynchronizationTag & tag) const { + AKANTU_DEBUG_IN(); + + UInt size = 0; + if (elements.empty()) { + return 0; + } + if (elements(0).kind() == _ek_regular) { + size = HeatTransferModel::getNbData(elements, tag); + } else if (elements(0).kind() == _ek_cohesive) { + UInt nb_nodes_per_element = 0; + Array::const_iterator it = elements.begin(); + Array::const_iterator end = elements.end(); + for (; it != end; ++it) { + const Element & el = *it; + nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); + } + switch (tag) { + case SynchronizationTag::_htm_temperature: { + size += nb_nodes_per_element * sizeof(Real); // temperature + break; + } + case SynchronizationTag::_htm_gradient_temperature: { + // temperature gradient + size += getNbIntegrationPoints(elements, "InterfacesFEEngine") * + spatial_dimension * sizeof(Real); + size += nb_nodes_per_element * sizeof(Real); // nodal temperatures + break; + } + default: { + AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); + } + } + } + AKANTU_DEBUG_OUT(); + return size; +} + +/* -------------------------------------------------------------------------- + */ +inline void +HeatTransferInterfaceModel::packData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) const { + if (elements.empty()) { + return; + } + if (elements(0).kind() == _ek_regular) { + HeatTransferModel::packData(buffer, elements, tag); + } else if (elements(0).kind() == _ek_cohesive) { + switch (tag) { + case SynchronizationTag::_htm_temperature: { + packNodalDataHelper(*temperature, buffer, elements, mesh); + break; + } + case SynchronizationTag::_htm_gradient_temperature: { + packElementalDataHelper(temperature_gradient, buffer, elements, true, + getFEEngine("InterfacesFEEngine")); + packNodalDataHelper(*temperature, buffer, elements, mesh); + break; + } + default: { + AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); + } + } + } +} + +/* -------------------------------------------------------------------------- + */ +inline void +HeatTransferInterfaceModel::unpackData(CommunicationBuffer & buffer, + const Array & elements, + const SynchronizationTag & tag) { + if (elements.empty()) { + return; + } + if (elements(0).kind() == _ek_regular) { + HeatTransferModel::unpackData(buffer, elements, tag); + } else if (elements(0).kind() == _ek_cohesive) { + switch (tag) { + case SynchronizationTag::_htm_temperature: { + unpackNodalDataHelper(*temperature, buffer, elements, mesh); + break; + } + case SynchronizationTag::_htm_gradient_temperature: { + unpackElementalDataHelper(temperature_gradient, buffer, elements, true, + getFEEngine("InterfacesFEEngine")); + unpackNodalDataHelper(*temperature, buffer, elements, mesh); + break; + } + default: { + AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); + } + } + } +} + +} // namespace akantu + +#endif /* AKANTU_HEAT_TRANSFER_MODEL_INLINE_IMPL_HH_ */