diff --git a/python/py_model.cc b/python/py_model.cc index 795b30ac9..f688183b9 100644 --- a/python/py_model.cc +++ b/python/py_model.cc @@ -1,89 +1,93 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_model(py::module & mod) { py::class_(mod, "DOFManager") .def("getMatrix", &DOFManager::getMatrix, py::return_value_policy::reference) .def( "getNewMatrix", [](DOFManager & self, const std::string & name, const std::string & matrix_to_copy_id) -> decltype(auto) { return self.getNewMatrix(name, matrix_to_copy_id); }, py::return_value_policy::reference) .def( "getResidual", [](DOFManager & self) -> decltype(auto) { return self.getResidual(); }, py::return_value_policy::reference) .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs) .def( - "hasMatrix", - [](DOFManager & self, const ID & name) -> bool { - return self.hasMatrix(name); - }, - py::arg("name")) + "hasMatrix", + [](DOFManager & self, const ID & name) -> bool { + return self.hasMatrix(name); + }, + py::arg("name")) .def("assembleToResidual", &DOFManager::assembleToResidual); py::class_(mod, "NonLinearSolver") .def( "set", [](NonLinearSolver & self, const std::string & id, const Real & val) { if (id == "max_iterations") { self.set(id, int(val)); } else { self.set(id, val); } }) .def("set", [](NonLinearSolver & self, const std::string & id, const SolveConvergenceCriteria & val) { self.set(id, val); }); py::class_(mod, "ModelSolver", py::multiple_inheritance()) .def("getNonLinearSolver", (NonLinearSolver & (ModelSolver::*)(const ID &)) & ModelSolver::getNonLinearSolver, py::arg("solver_id") = "", py::return_value_policy::reference) .def("solveStep", [](ModelSolver & self) { self.solveStep(); }) .def("solveStep", [](ModelSolver & self, const ID & solver_id) { self.solveStep(solver_id); }); py::class_(mod, "Model", py::multiple_inheritance()) .def("setBaseName", &Model::setBaseName) .def("getFEEngine", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("getFEEngineBoundary", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("addDumpFieldVector", &Model::addDumpFieldVector) .def("addDumpField", &Model::addDumpField) .def("setBaseNameToDumper", &Model::setBaseNameToDumper) .def("addDumpFieldVectorToDumper", &Model::addDumpFieldVectorToDumper) .def("addDumpFieldToDumper", &Model::addDumpFieldToDumper) .def("dump", py::overload_cast<>(&Model::dump)) .def("dump", py::overload_cast(&Model::dump)) .def("dump", py::overload_cast(&Model::dump)) + .def("dump", py::overload_cast(&Model::dump)) + .def("dump", py::overload_cast(&Model::dump)) + .def("dump", + py::overload_cast(&Model::dump)) .def("initNewSolver", &Model::initNewSolver) .def("getDOFManager", &Model::getDOFManager, py::return_value_policy::reference) .def("assembleMatrix", &Model::assembleMatrix); } } // namespace akantu diff --git a/python/py_solid_mechanics_model.cc b/python/py_solid_mechanics_model.cc index b92975469..e5f6d28b2 100644 --- a/python/py_solid_mechanics_model.cc +++ b/python/py_solid_mechanics_model.cc @@ -1,117 +1,110 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ void register_solid_mechanics_model(py::module & mod) { py::class_(mod, "SolidMechanicsModelOptions") .def(py::init(), py::arg("_analysis_method") = _explicit_lumped_mass); py::class_(mod, "SolidMechanicsModel", py::multiple_inheritance()) .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "solid_mechanics_model", py::arg("model_type") = ModelType::_solid_mechanics_model) .def( "initFull", [](SolidMechanicsModel & self, const SolidMechanicsModelOptions & options) { self.initFull(options); }, py::arg("option") = SolidMechanicsModelOptions()) .def( "initFull", [](SolidMechanicsModel & self, const AnalysisMethod & analysis_method) { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method")) .def_deprecated("applyDirichletBC", "Deprecated: use applyBC") .def("applyBC", [](SolidMechanicsModel & self, BC::Dirichlet::DirichletFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("applyBC", [](SolidMechanicsModel & self, BC::Neumann::NeumannFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("setTimeStep", &SolidMechanicsModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def("getEnergy", py::overload_cast( &SolidMechanicsModel::getEnergy), py::arg("energy_id")) .def("getEnergy", py::overload_cast( &SolidMechanicsModel::getEnergy), py::arg("energy_id"), py::arg("group_id")) .def_function(assembleStiffnessMatrix) .def_function(assembleInternalForces) .def_function(assembleMass) .def_function(assembleMassLumped) .def_function(getStableTimeStep) .def_function_nocopy(getExternalForce) .def_function_nocopy(getDisplacement) .def_function_nocopy(getPreviousDisplacement) .def_function_nocopy(getIncrement) .def_function_nocopy(getInternalForce) .def_function_nocopy(getMass) .def_function_nocopy(getVelocity) .def_function_nocopy(getAcceleration) .def_function_nocopy(getInternalForce) .def_function_nocopy(getBlockedDOFs) .def_function_nocopy(getMesh) - .def("dump", py::overload_cast<>(&SolidMechanicsModel::dump)) - .def("dump", - py::overload_cast(&SolidMechanicsModel::dump)) - .def("dump", py::overload_cast( - &SolidMechanicsModel::dump)) - .def("dump", py::overload_cast( - &SolidMechanicsModel::dump)) .def("getMaterial", py::overload_cast(&SolidMechanicsModel::getMaterial), py::return_value_policy::reference) .def("getMaterial", py::overload_cast( &SolidMechanicsModel::getMaterial), py::return_value_policy::reference) .def("getMaterialIndex", &SolidMechanicsModel::getMaterialIndex) .def("setMaterialSelector", &SolidMechanicsModel::setMaterialSelector) .def("getMaterialSelector", &SolidMechanicsModel::getMaterialSelector); } } // namespace akantu diff --git a/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index bbb161bf5..1035d1907 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,945 +1,918 @@ /** * @file heat_transfer_model.cc * * @author Guillaume Anciaux * @author Lucas Frerot * @author Emil Gallyamov * @author David Simon Kammer * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Tue Feb 20 2018 * * @brief Implementation of HeatTransferModel class * * * Copyright (©) 2010-2018 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 "dumpable_inline_impl.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_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { namespace heat_transfer { namespace details { class ComputeRhoFunctor { public: ComputeRhoFunctor(const HeatTransferModel & model) : model(model){}; void operator()(Matrix & rho, const Element & /*unused*/) { rho.set(model.getCapacity() * model.getDensity()); } private: const HeatTransferModel & model; }; } // namespace details } // namespace heat_transfer /* -------------------------------------------------------------------------- */ HeatTransferModel::HeatTransferModel(Mesh & mesh, UInt dim, const ID & id) : Model(mesh, ModelType::_heat_transfer_model, dim, id), 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) { AKANTU_DEBUG_IN(); conductivity = Matrix(this->spatial_dimension, this->spatial_dimension); this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_htm_temperature); this->registerSynchronizer(synchronizer, SynchronizationTag::_htm_gradient_temperature); } registerFEEngineObject(id + ":fem", mesh, spatial_dimension); #ifdef AKANTU_USE_IOHELPER 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() { auto & fem = this->getFEEngine(); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); temperature_on_qpoints.initialize(fem, _nb_component = 1); temperature_gradient.initialize(fem, _nb_component = spatial_dimension); conductivity_on_qpoints.initialize(fem, _nb_component = spatial_dimension * spatial_dimension); k_gradt_on_qpoints.initialize(fem, _nb_component = spatial_dimension); } /* -------------------------------------------------------------------------- */ FEEngine & HeatTransferModel::getFEEngineBoundary(const ID & name) { return aka::as_type(getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ HeatTransferModel::~HeatTransferModel() = default; /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor compute_rho(*this); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldLumped(compute_rho, "M", "temperature", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType HeatTransferModel::getMatrixType(const ID & matrix_id) { if (matrix_id == "K" or matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleConductivityMatrix(); } else if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacity(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacityLumped(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleResidual() { AKANTU_DEBUG_IN(); this->assembleInternalHeatRate(); this->getDOFManager().assembleToResidual("temperature", *this->external_heat_rate, 1); this->getDOFManager().assembleToResidual("temperature", *this->internal_heat_rate, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::predictor() { ++temperature_release; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped() { AKANTU_DEBUG_IN(); if (!this->getDOFManager().hasLumpedMatrix("M")) { this->getDOFManager().getNewLumpedMatrix("M"); } this->getDOFManager().zeroLumpedMatrix("M"); assembleCapacityLumped(_not_ghost); assembleCapacityLumped(_ghost); need_to_reassemble_capacity_lumped = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->temperature, 1, "temperature"); this->allocNodalField(this->external_heat_rate, 1, "external_heat_rate"); this->allocNodalField(this->internal_heat_rate, 1, "internal_heat_rate"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); if (!dof_manager.hasDOFs("temperature")) { dof_manager.registerDOFs("temperature", *this->temperature, _dst_nodal); dof_manager.registerBlockedDOFs("temperature", *this->blocked_dofs); } if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(this->temperature_rate, 1, "temperature_rate"); if (!dof_manager.hasDOFsDerivatives("temperature", 1)) { dof_manager.registerDOFsDerivative("temperature", 1, *this->temperature_rate); } } } /* -------------------------------------------------------------------------- */ std::tuple HeatTransferModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions HeatTransferModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_forward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_pseudo_time; options.solution_type["temperature"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_forward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; } else { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_backward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleConductivityMatrix() { AKANTU_DEBUG_IN(); this->computeConductivityOnQuadPoints(_not_ghost); if (conductivity_release[_not_ghost] == conductivity_matrix_release) { return; } AKANTU_DEBUG_ASSERT(this->getDOFManager().hasMatrix("K"), "The K matrix has not been initialized yet."); this->getDOFManager().zeroMatrix("K"); auto & fem = this->getFEEngine(); for (auto && type : mesh.elementTypes(spatial_dimension)) { auto nb_element = mesh.getNbElement(type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); fem.computeBtDB(conductivity_on_qpoints(type), *bt_d_b, 2, type); /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto K_e = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_e"); fem.integrate(*bt_d_b, *K_e, nb_nodes_per_element * nb_nodes_per_element, type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "temperature", *K_e, type, _not_ghost, _symmetric); } conductivity_matrix_release = conductivity_release[_not_ghost]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeConductivityOnQuadPoints(GhostType ghost_type) { // if already computed once check if need to compute if (not initial_conductivity[ghost_type]) { // if temperature did not change, conductivity will not vary if (temperature_release == conductivity_release[ghost_type]) { return; } // if conductivity_variation is 0 no need to recompute if (conductivity_variation == 0.) { return; } } for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & temperature_interpolated = temperature_on_qpoints(type, ghost_type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type, ghost_type); auto & cond = conductivity_on_qpoints(type, ghost_type); for (auto && tuple : zip(make_view(cond, spatial_dimension, spatial_dimension), temperature_interpolated)) { auto & C = std::get<0>(tuple); auto & T = std::get<1>(tuple); C = conductivity; Matrix variation(spatial_dimension, spatial_dimension, conductivity_variation * (T - T_ref)); // @TODO: Guillaume are you sure ? why due you compute variation then ? C += conductivity_variation; } } conductivity_release[ghost_type] = temperature_release; initial_conductivity[ghost_type] = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeKgradT(GhostType ghost_type) { computeConductivityOnQuadPoints(ghost_type); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & gradient = temperature_gradient(type, ghost_type); this->getFEEngine().gradientOnIntegrationPoints(*temperature, gradient, 1, type, ghost_type); for (auto && values : zip(make_view(conductivity_on_qpoints(type, ghost_type), spatial_dimension, spatial_dimension), make_view(gradient, spatial_dimension), make_view(k_gradt_on_qpoints(type, ghost_type), spatial_dimension))) { const auto & C = std::get<0>(values); const auto & BT = std::get<1>(values); auto & k_BT = std::get<2>(values); k_BT.mul(C, BT); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleInternalHeatRate() { AKANTU_DEBUG_IN(); this->internal_heat_rate->zero(); this->synchronize(SynchronizationTag::_htm_temperature); auto & fem = this->getFEEngine(); for (auto ghost_type : ghost_types) { // compute k \grad T computeKgradT(ghost_type); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto & k_gradt_on_qpoints_vect = k_gradt_on_qpoints(type, ghost_type); UInt nb_quad_points = k_gradt_on_qpoints_vect.size(); Array bt_k_gT(nb_quad_points, nb_nodes_per_element); fem.computeBtD(k_gradt_on_qpoints_vect, bt_k_gT, type, ghost_type); UInt nb_elements = mesh.getNbElement(type, ghost_type); Array int_bt_k_gT(nb_elements, nb_nodes_per_element); fem.integrate(bt_k_gT, int_bt_k_gT, nb_nodes_per_element, type, ghost_type); this->getDOFManager().assembleElementalArrayLocalArray( int_bt_k_gT, *this->internal_heat_rate, type, ghost_type, -1); } } 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); } } for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array coord(0, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, _not_ghost); auto el_coord = coord.begin(spatial_dimension, nb_nodes_per_element); UInt nb_element = mesh.getNbElement(type); for (UInt el = 0; el < nb_element; ++el, ++el_coord) { el_size = getFEEngine().getElementInradius(*el_coord, type); 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 / 4. * density * capacity / conductivitymax; mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("heat_transfer").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void HeatTransferModel::readMaterials() { auto sect = this->getParserSection(); if (not std::get<1>(sect)) { const auto & section = std::get<0>(sect); this->parseSection(section); } conductivity_on_qpoints.set(conductivity); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); readMaterials(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacity() { AKANTU_DEBUG_IN(); auto ghost_type = _not_ghost; this->getDOFManager().zeroMatrix("M"); auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor rho_functor(*this); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldMatrix(rho_functor, "M", "temperature", this->getDOFManager(), type, ghost_type); } need_to_reassemble_capacity = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeRho(Array & rho, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); FEEngine & fem = this->getFEEngine(); UInt nb_element = mesh.getNbElement(type, ghost_type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); rho.resize(nb_element * nb_quadrature_points); rho.set(this->capacity); // 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(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::computeThermalEnergyByNode() { AKANTU_DEBUG_IN(); Real ethermal = 0.; for (auto && pair : enumerate(make_view( *internal_heat_rate, internal_heat_rate->getNbComponent()))) { auto n = std::get<0>(pair); auto & heat_rate = std::get<1>(pair); Real heat = 0.; bool is_local_node = mesh.isLocalOrMasterNode(n); bool count_node = is_local_node; for (UInt i = 0; i < heat_rate.size(); ++i) { if (count_node) { heat += heat_rate[i] * time_step; } } ethermal += heat; } mesh.getCommunicator().allReduce(ethermal, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ethermal; } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::getThermalEnergy( iterator Eth, Array::const_iterator T_it, const Array::const_iterator & T_end) const { for (; T_it != T_end; ++T_it, ++Eth) { *Eth = capacity * density * *T_it; } } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy(ElementType type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Vector Eth_on_quarature_points(nb_quadrature_points); auto T_it = this->temperature_on_qpoints(type).begin(); T_it += index * nb_quadrature_points; auto 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; auto & fem = getFEEngine(); for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_element = mesh.getNbElement(type, _not_ghost); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, _not_ghost); Array Eth_per_quad(nb_element * nb_quadrature_points, 1); auto & temperature_interpolated = temperature_on_qpoints(type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type); auto T_it = temperature_interpolated.begin(); auto T_end = temperature_interpolated.end(); getThermalEnergy(Eth_per_quad.begin(), T_it, T_end); Eth += fem.integrate(Eth_per_quad, type); } 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 mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id, 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 std::shared_ptr 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.get(); auto field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr HeatTransferModel::createNodalFieldReal( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { if (field_name == "capacity_lumped") { AKANTU_EXCEPTION( "Capacity lumped is a nodal field now stored in the DOF manager." "Therefore it cannot be used by a dumper anymore"); } std::map *> real_nodal_fields; real_nodal_fields["temperature"] = temperature.get(); real_nodal_fields["temperature_rate"] = temperature_rate.get(); real_nodal_fields["external_heat_rate"] = external_heat_rate.get(); real_nodal_fields["internal_heat_rate"] = internal_heat_rate.get(); real_nodal_fields["increment"] = increment.get(); std::shared_ptr field = mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr HeatTransferModel::createElementalField( const std::string & field_name, const std::string & group_name, bool /*padding_flag*/, UInt /*spatial_dimension*/, ElementKind element_kind) { std::shared_ptr field; 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); 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); field = mesh.createElementalField( conductivity_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ std::shared_ptr HeatTransferModel::createElementalField( - __attribute__((unused)) const std::string & field_name, - __attribute__((unused)) const std::string & group_name, - __attribute__((unused)) bool padding_flag, - __attribute__((unused)) ElementKind element_kind) { + const std::string & /* field_name*/, const std::string & /*group_name*/, + bool /*padding_flag*/, ElementKind /*element_kind*/) { return nullptr; } /* -------------------------------------------------------------------------- */ -std::shared_ptr HeatTransferModel::createNodalFieldBool( - __attribute__((unused)) const std::string & field_name, - __attribute__((unused)) const std::string & group_name, - __attribute__((unused)) bool padding_flag) { +std::shared_ptr +HeatTransferModel::createNodalFieldBool(const std::string & /*field_name*/, + const std::string & /*group_name*/, + bool /*padding_flag*/) { return nullptr; } /* -------------------------------------------------------------------------- */ -std::shared_ptr HeatTransferModel::createNodalFieldReal( - __attribute__((unused)) const std::string & field_name, - __attribute__((unused)) const std::string & group_name, - __attribute__((unused)) bool padding_flag) { +std::shared_ptr +HeatTransferModel::createNodalFieldReal(const std::string & /*field_name*/, + const std::string & /*group_name*/, + bool /*padding_flag*/) { return nullptr; } #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); } - /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbData(const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = indexes.size(); switch (tag) { case SynchronizationTag::_htm_temperature: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case SynchronizationTag::_htm_temperature: { buffer << (*temperature)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case SynchronizationTag::_htm_temperature: { buffer >> (*temperature)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::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; 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) * 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 HeatTransferModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { switch (tag) { case SynchronizationTag::_htm_temperature: { packNodalDataHelper(*temperature, buffer, elements, mesh); break; } case SynchronizationTag::_htm_gradient_temperature: { packElementalDataHelper(temperature_gradient, buffer, elements, true, getFEEngine()); packNodalDataHelper(*temperature, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { switch (tag) { case SynchronizationTag::_htm_temperature: { unpackNodalDataHelper(*temperature, buffer, elements, mesh); break; } case SynchronizationTag::_htm_gradient_temperature: { unpackElementalDataHelper(temperature_gradient, buffer, elements, true, getFEEngine()); unpackNodalDataHelper(*temperature, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/heat_transfer/heat_transfer_model.hh b/src/model/heat_transfer/heat_transfer_model.hh index aca1ad271..0ba0ef303 100644 --- a/src/model/heat_transfer/heat_transfer_model.hh +++ b/src/model/heat_transfer/heat_transfer_model.hh @@ -1,327 +1,315 @@ /** * @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: Mon Feb 05 2018 * * @brief Model of Heat Transfer * * * Copyright (©) 2010-2018 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 "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_HEAT_TRANSFER_MODEL_HH_ #define AKANTU_HEAT_TRANSFER_MODEL_HH_ namespace akantu { template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu namespace akantu { class HeatTransferModel : public Model, public DataAccessor, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using FEEngineType = FEEngineTemplate; HeatTransferModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "heat_transfer_model"); ~HeatTransferModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// generic function to initialize everything ready for explicit dynamics void initFullImpl(const ModelOptions & options) override; /// read one material file to instantiate all the materials void readMaterials(); /// allocate all vectors void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; /// initialize the model void initModel() override; void predictor() override; /// compute the heat flux void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback to assemble a Matrix void assembleMatrix(const ID & matrix_id) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override; std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; /* ------------------------------------------------------------------------ */ /* Methods for explicit */ /* ------------------------------------------------------------------------ */ public: /// compute and get the stable time step Real getStableTimeStep(); /// set the stable timestep void setTimeStep(Real time_step, const ID & solver_id = "") override; // temporary protection to prevent bad usage: should check for bug protected: /// compute the internal heat flux \todo Need code review: currently not /// public method void assembleInternalHeatRate(); public: /// calculate the lumped capacity vector for heat transfer problem void assembleCapacityLumped(); public: /// assemble the conductivity matrix void assembleConductivityMatrix(); /// assemble the conductivity matrix void assembleCapacity(); /// compute the capacity on quadrature points void computeRho(Array & rho, ElementType type, GhostType ghost_type); private: /// calculate the lumped capacity vector for heat transfer problem (w /// ghost type) void assembleCapacityLumped(GhostType ghost_type); /// compute the conductivity tensor for each quadrature point in an array void computeConductivityOnQuadPoints(GhostType ghost_type); /// compute vector \f[k \grad T\f] for each quadrature point void computeKgradT(GhostType ghost_type); /// compute the thermal energy Real computeThermalEnergyByNode(); /* ------------------------------------------------------------------------ */ /* 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; inline UInt getNbData(const Array & indexes, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; - 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); - - void dump() override; - - virtual void dump(UInt step); - - virtual void dump(Real time, UInt step); - /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: 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); /// get the assembled heat flux AKANTU_GET_MACRO(InternalHeatRate, *internal_heat_rate, Array &); /// get the boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the external heat rate vector AKANTU_GET_MACRO(ExternalHeatRate, *external_heat_rate, Array &); /// get the temperature gradient 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); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureOnQpoints, temperature_on_qpoints, Real); /// internal variables AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(KgradT, k_gradt_on_qpoints, Real); /// get the temperature AKANTU_GET_MACRO(Temperature, *temperature, Array &); /// get the temperature derivative AKANTU_GET_MACRO(TemperatureRate, *temperature_rate, Array &); /// get the energy denominated by thermal Real getEnergy(const std::string & energy_id, 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(ElementType type, UInt index); /// get the thermal energy for a given element Real getThermalEnergy(); protected: /* ------------------------------------------------------------------------ */ FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ----------------------------------------------------------------------- */ template void getThermalEnergy(iterator Eth, Array::const_iterator T_it, const Array::const_iterator & T_end) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// temperatures array std::unique_ptr> temperature; /// temperatures derivatives array std::unique_ptr> temperature_rate; /// increment array (@f$\delta \dot T@f$ or @f$\delta T@f$) std::unique_ptr> increment; /// 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 \f[k \grad T\f] on quad points ElementTypeMapArray k_gradt_on_qpoints; /// external flux vector std::unique_ptr> external_heat_rate; /// residuals array std::unique_ptr> internal_heat_rate; /// boundary vector std::unique_ptr> blocked_dofs; // realtime // Real time; /// capacity Real capacity; // conductivity matrix Matrix 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 // Real conductivitymax; bool need_to_reassemble_capacity{true}; bool need_to_reassemble_capacity_lumped{true}; UInt temperature_release{0}; UInt conductivity_matrix_release{UInt(-1)}; std::unordered_map initial_conductivity{{_not_ghost, true}, {_ghost, true}}; std::unordered_map conductivity_release{{_not_ghost, 0}, {_ghost, 0}}; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_model_inline_impl.hh" #endif /* AKANTU_HEAT_TRANSFER_MODEL_HH_ */ diff --git a/src/model/model.cc b/src/model/model.cc index cf8b86a54..e1180198d 100644 --- a/src/model/model.cc +++ b/src/model/model.cc @@ -1,339 +1,364 @@ /** * @file model.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Mon Oct 03 2011 * @date last modification: Tue Feb 20 2018 * * @brief implementation of model common parts * * * Copyright (©) 2010-2018 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 "model.hh" #include "communicator.hh" #include "data_accessor.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Model::Model(Mesh & mesh, const ModelType & type, UInt dim, const ID & id) : ModelSolver(mesh, type, id), mesh(mesh), spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension() : dim), parser(getStaticParser()) { AKANTU_DEBUG_IN(); this->mesh.registerEventHandler(*this, _ehp_model); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Model::~Model() = default; /* -------------------------------------------------------------------------- */ void Model::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); method = options.analysis_method; if (!this->hasDefaultSolver()) { this->initNewSolver(this->method); } initModel(); initFEEngineBoundary(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Model::initNewSolver(const AnalysisMethod & method) { ID solver_name; TimeStepSolverType tss_type; std::tie(solver_name, tss_type) = this->getDefaultSolverID(method); if (not this->hasSolver(solver_name)) { ModelSolverOptions options = this->getDefaultSolverOptions(tss_type); this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type); for (auto && is_type : options.integration_scheme_type) { if (!this->hasIntegrationScheme(solver_name, is_type.first)) { this->setIntegrationScheme(solver_name, is_type.first, is_type.second, options.solution_type[is_type.first]); } } } this->method = method; this->setDefaultSolver(solver_name); } /* -------------------------------------------------------------------------- */ void Model::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name, const std::string & dumper_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup() { for (auto & group : mesh.iterateElementGroups()) { group.dump(); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory) { for (auto & group : mesh.iterateElementGroups()) { group.setDirectory(directory); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Model::setGroupBaseName(const std::string & basename, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setBaseName(basename); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Model::getGroupDumper(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); return group.getDumper(); } /* -------------------------------------------------------------------------- */ // DUMPER stuff /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & field_id, std::shared_ptr field, DumperIOHelper & dumper) { #ifdef AKANTU_USE_IOHELPER dumper.registerField(field_id, std::move(field)); #endif } /* -------------------------------------------------------------------------- */ void Model::addDumpField(const std::string & field_id) { this->addDumpFieldToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldVector(const std::string & field_id) { this->addDumpFieldVectorToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldTensor(const std::string & field_id) { this->addDumpFieldTensorToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::setBaseName(const std::string & field_id) { mesh.setBaseName(field_id); } /* -------------------------------------------------------------------------- */ void Model::setBaseNameToDumper(const std::string & dumper_name, const std::string & basename) { mesh.setBaseNameToDumper(dumper_name, basename); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", dumper_default_element_kind, false); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id, group_name, dumper_default_element_kind, false); } /* -------------------------------------------------------------------------- */ void Model::removeDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void Model::removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.removeDumpFieldFromDumper(dumper_name, field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", dumper_default_element_kind, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name, dumper_default_element_kind, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", dumper_default_element_kind, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag) { this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name, this->spatial_dimension, element_kind, padding_flag); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, ElementKind element_kind, bool padding_flag) { #ifdef AKANTU_USE_IOHELPER std::shared_ptr field; if (!field) { field = this->createNodalFieldReal(field_id, group_name, padding_flag); } if (!field) { field = this->createNodalFieldUInt(field_id, group_name, padding_flag); } if (!field) { field = this->createNodalFieldBool(field_id, group_name, padding_flag); } if (!field) { field = this->createElementalField(field_id, group_name, padding_flag, spatial_dimension, element_kind); } if (!field) { field = this->mesh.createFieldFromAttachedData(field_id, group_name, element_kind); } if (!field) { field = this->mesh.createFieldFromAttachedData(field_id, group_name, element_kind); } #ifndef AKANTU_NDEBUG if (!field) { AKANTU_DEBUG_WARNING("No field could be found based on name: " << field_id); } #endif if (field) { DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name, group_name); this->addDumpGroupFieldToDumper(field_id, field, dumper); } #endif } /* -------------------------------------------------------------------------- */ -void Model::dump() { mesh.dump(); } +void Model::dump(const std::string & dumper_name) { + mesh.dump(dumper_name); +} /* -------------------------------------------------------------------------- */ -void Model::dump(UInt step) { mesh.dump(step); } +void Model::dump(const std::string & dumper_name, UInt step) { + mesh.dump(dumper_name, step); +} + +/* ------------------------------------------------------------------------- */ +void Model::dump(const std::string & dumper_name, Real time, + UInt step) { + mesh.dump(dumper_name, time, step); +} /* -------------------------------------------------------------------------- */ -void Model::dump(Real time, UInt step) { mesh.dump(time, step); } +void Model::dump() { + auto default_dumper = mesh.getDefaultDumperName(); + this->dump(default_dumper); +} + +/* -------------------------------------------------------------------------- */ +void Model::dump(UInt step) { + auto default_dumper = mesh.getDefaultDumperName(); + this->dump(default_dumper, step); +} + +/* -------------------------------------------------------------------------- */ +void Model::dump(Real time, UInt step) { + auto default_dumper = mesh.getDefaultDumperName(); + this->dump(default_dumper, time, step); +} /* -------------------------------------------------------------------------- */ void Model::setDirectory(const std::string & directory) { mesh.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Model::setDirectoryToDumper(const std::string & dumper_name, const std::string & directory) { mesh.setDirectoryToDumper(dumper_name, directory); } /* -------------------------------------------------------------------------- */ void Model::setTextModeToDumper() { mesh.setTextModeToDumper(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/model.hh b/src/model/model.hh index a46846081..c23ca4807 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,351 +1,357 @@ /** * @file model.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief Interface of a model * * * Copyright (©) 2010-2018 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 "aka_named_argument.hh" #include "fe_engine.hh" #include "mesh.hh" #include "model_options.hh" #include "model_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MODEL_HH_ #define AKANTU_MODEL_HH_ namespace akantu { class SynchronizerRegistry; class Parser; class DumperIOHelper; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { class Model : public ModelSolver, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Model(Mesh & mesh, const ModelType & type, UInt dim = _all_dimensions, const ID & id = "model"); ~Model() override; using FEEngineMap = std::map>; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: virtual void initFullImpl(const ModelOptions & options); public: template std::enable_if_t::value> initFull(pack &&... _pack) { switch (this->model_type) { #ifdef AKANTU_SOLID_MECHANICS case ModelType::_solid_mechanics_model: this->initFullImpl(SolidMechanicsModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_COHESIVE_ELEMENT case ModelType::_solid_mechanics_model_cohesive: this->initFullImpl(SolidMechanicsModelCohesiveOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_HEAT_TRANSFER case ModelType::_heat_transfer_model: this->initFullImpl(HeatTransferModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_EMBEDDED case ModelType::_embedded_model: this->initFullImpl(EmbeddedInterfaceModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif default: this->initFullImpl(ModelOptions{use_named_args, std::forward(_pack)...}); } } template std::enable_if_t::value> initFull(pack &&... _pack) { this->initFullImpl(std::forward(_pack)...); } /// initialize a new solver if needed void initNewSolver(const AnalysisMethod & method); protected: /// get some default values for derived classes virtual std::tuple getDefaultSolverID(const AnalysisMethod & method) = 0; virtual void initModel() = 0; virtual void initFEEngineBoundary(); /// function to print the containt of the class void printself(std::ostream & /*stream*/, int /*indent*/ = 0) const override{}; public: /* ------------------------------------------------------------------------ */ /* Access to the dumpable interface of the boundaries */ /* ------------------------------------------------------------------------ */ /// Dump the data for a given group void dumpGroup(const std::string & group_name); void dumpGroup(const std::string & group_name, const std::string & dumper_name); /// Dump the data for all boundaries void dumpGroup(); /// Set the directory for a given group void setGroupDirectory(const std::string & directory, const std::string & group_name); /// Set the directory for all boundaries void setGroupDirectory(const std::string & directory); /// Set the base name for a given group void setGroupBaseName(const std::string & basename, const std::string & group_name); /// Get the internal dumper of a given group DumperIOHelper & getGroupDumper(const std::string & group_name); /* ------------------------------------------------------------------------ */ /* Function for non local capabilities */ /* ------------------------------------------------------------------------ */ virtual void updateDataForNonLocalCriterion(__attribute__((unused)) ElementTypeMapReal & criterion) { AKANTU_TO_IMPLEMENT(); } protected: template void allocNodalField(std::unique_ptr> & array, UInt nb_component, const ID & name) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get id of model AKANTU_GET_MACRO(ID, id, const ID &) /// get the number of surfaces AKANTU_GET_MACRO(Mesh, mesh, Mesh &) /// synchronize the boundary in case of parallel run virtual void synchronizeBoundaries(){}; /// return the fem object associated with a provided name inline FEEngine & getFEEngine(const ID & name = "") const; /// return the fem boundary object associated with a provided name virtual FEEngine & getFEEngineBoundary(const ID & name = ""); /// register a fem object associated with name template inline void registerFEEngineObject(const std::string & name, Mesh & mesh, UInt spatial_dimension); /// unregister a fem object associated with name inline void unRegisterFEEngineObject(const std::string & name); /// return the synchronizer registry SynchronizerRegistry & getSynchronizerRegistry(); /// return the fem object associated with a provided name template inline FEEngineClass & getFEEngineClass(std::string name = "") const; /// return the fem boundary object associated with a provided name template inline FEEngineClass & getFEEngineClassBoundary(std::string name = ""); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /* ------------------------------------------------------------------------ */ - /* Pack and unpack helper functions */ + /* Pack and unpack hexlper functions */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbIntegrationPoints(const Array & elements, const ID & fem_id = ID()) const; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ void setTextModeToDumper(); virtual void addDumpGroupFieldToDumper(const std::string & field_id, std::shared_ptr field, DumperIOHelper & dumper); virtual void addDumpField(const std::string & field_id); virtual void addDumpFieldVector(const std::string & field_id); virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensor(const std::string & field_id); virtual void setBaseName(const std::string & field_id); virtual void setBaseNameToDumper(const std::string & dumper_name, const std::string & basename); virtual void addDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, ElementKind element_kind, bool padding_flag); virtual void removeDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual std::shared_ptr createNodalFieldReal(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createNodalFieldUInt(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createNodalFieldBool(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } virtual std::shared_ptr createElementalField(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/, UInt /*spatial_dimension*/, ElementKind /*kind*/) { return nullptr; } void setDirectory(const std::string & directory); void setDirectoryToDumper(const std::string & dumper_name, const std::string & directory); + + /* ------------------------------------------------------------------------ */ + 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 void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: friend std::ostream & operator<<(std::ostream & /*stream*/, const Model & /*_this*/); ID id; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// Mesh Mesh & mesh; /// Spatial dimension of the problem UInt spatial_dimension; /// the main fem object present in all models FEEngineMap fems; /// the fem object present in all models for boundaries FEEngineMap fems_boundary; /// default fem object std::string default_fem; /// parser to the pointer to use Parser & parser; /// default ElementKind for dumper ElementKind dumper_default_element_kind{_ek_regular}; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Model & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "model_inline_impl.hh" #endif /* AKANTU_MODEL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index d9c1f019d..48714facf 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,580 +1,580 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * * Copyright (©) 2010-2018 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 "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager_callback.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_HH_ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "solid_mechanics_model", ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; public: /// initialize all internal arrays for materials virtual void initMaterials(); protected: /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(ElementType type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(GhostType ghost_type) override; void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, ElementKind element_kind); //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, ElementKind kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, ElementKind kind, GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(ElementKind kind); std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; - 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); + void dump(const std::string & dumper_name) override; + void dump(const std::string & dumper_name, UInt step) override; + void dump(const std::string & dumper_name, Real time, UInt step) override; void dump() override; void dump(UInt step) override; void dump(Real time, UInt step) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement); /// get the SolidMechanicsModel::previous_displacement array AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement); /// get the SolidMechanicsModel::current_position array const Array & getCurrentPosition(); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment); /// get the lumped SolidMechanicsModel::mass array AKANTU_GET_MACRO_DEREF_PTR(Mass, mass); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force); /// get the SolidMechanicsModel::force array (external forces) [[deprecated("Use getExternalForce instead of this function")]] Array & getForce() { return getExternalForce(); } /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force); /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs); /// get an iterable on the materials inline decltype(auto) getMaterials(); /// get an iterable on the materials inline decltype(auto) getMaterials() const; /// get a particular material (by numerical material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by numerical material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for an element Real getEnergy(const std::string & energy_id, ElementType type, UInt index); /// compute the energy for an element Real getEnergy(const std::string & energy_id, const Element & element) { return getEnergy(energy_id, element.type, element.element); } /// compute the energy for an element group Real getEnergy(const ID & energy_id, const ID & group_id); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, // material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector, MaterialSelector &); void setMaterialSelector(std::shared_ptr material_selector) { this->material_selector = std::move(material_selector); } /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: /// compute the stable time step Real getStableTimeStep(GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// release version of the displacement array UInt displacement_release{0}; /// release version of the current_position array UInt current_position_release{0}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// mapping between material name and material internal id std::map materials_names_to_id; protected: /// conversion coefficient form force/mass to acceleration Real f_m2a{1.0}; /// displacements array std::unique_ptr> displacement; /// displacements array at the previous time step (used in finite deformation) std::unique_ptr> previous_displacement; /// increment of displacement std::unique_ptr> displacement_increment; /// lumped mass array std::unique_ptr> mass; /// velocities array std::unique_ptr> velocity; /// accelerations array std::unique_ptr> acceleration; /// external forces array std::unique_ptr> external_force; /// internal forces array std::unique_ptr> internal_force; /// array specifing if a degree of freedom is blocked or not std::unique_ptr> blocked_dofs; /// array of current position used during update residual std::unique_ptr> current_position; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; /// list of used materials std::vector> materials; /// class defining of to choose a material std::shared_ptr material_selector; using flatten_internal_map = std::map, std::unique_ptr>>; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; /// tells if the material are instantiated bool are_materials_instantiated{false}; friend class Material; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.hh" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */ diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh index 67a05709f..fd2ae28bc 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.hh +++ b/src/model/structural_mechanics/structural_mechanics_model.hh @@ -1,338 +1,338 @@ /** * @file structural_mechanics_model.hh * * @author Fabian Barras * @author Lucas Frerot * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Tue Feb 20 2018 * * @brief Particular implementation of the structural elements in the * StructuralMechanicsModel * * * Copyright (©) 2010-2018 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_named_argument.hh" #include "boundary_condition.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ #define AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeStructural; } // namespace akantu namespace akantu { struct StructuralMaterial { Real E{0}; Real A{1}; Real I{0}; Real Iz{0}; Real Iy{0}; Real GJ{0}; Real rho{0}; Real t{0}; Real nu{0}; }; class StructuralMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineType = FEEngineTemplate; StructuralMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "structural_mechanics_model"); ~StructuralMechanicsModel() override; /// Init full model void initFullImpl(const ModelOptions & options) override; /// Init boundary FEEngine void initFEEngineBoundary() override; /* ------------------------------------------------------------------------ */ /* Virtual methods from SolverCallback */ /* ------------------------------------------------------------------------ */ /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback to assemble a Matrix void assembleMatrix(const ID & matrix_id) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback to assemble the residual (rhs) void assembleResidual() override; void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return false; } /// compute kinetic energy Real getKineticEnergy(); /// compute potential energy Real getPotentialEnergy(); /// compute the specified energy Real getEnergy(const ID & energy); /* ------------------------------------------------------------------------ */ /* Virtual methods from Model */ /* ------------------------------------------------------------------------ */ protected: /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; static UInt getNbDegreeOfFreedom(ElementType type); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; /// initialize the model void initModel() override; /// compute the stresses per elements void computeStresses(); /// compute the nodal forces void assembleInternalForce(); /// compute the nodal forces for an element type void assembleInternalForce(ElementType type, GhostType gt); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// assemble the mass matrix for consistent mass resolutions void assembleMassMatrix(); protected: /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMassMatrix(GhostType ghost_type); /// computes rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// finish the computation of residual to solve in increment void updateResidualInternal(); /* ------------------------------------------------------------------------ */ private: template void assembleStiffnessMatrix(); template void computeStressOnQuad(); template void computeTangentModuli(Array & tangent_moduli); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; 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: /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the StructuralMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement_rotation, Array &); /// get the StructuralMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the StructuralMechanicsModel::acceleration vector, updated /// by /// StructuralMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the StructuralMechanicsModel::external_force vector AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the StructuralMechanicsModel::internal_force vector (boundary forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the StructuralMechanicsModel::boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(RotationMatrix, rotation_matrix, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Set_ID, set_ID, UInt); /** * \brief This function adds the `StructuralMaterial` material to the list of * materials managed by *this. * * It is important that this function might invalidate all references to * structural materials, that were previously optained by `getMaterial()`. * * \param material The new material. * * \return The ID of the material that was added. * * \note The return type is is new. */ UInt addMaterial(StructuralMaterial & material, const ID & name = ""); const StructuralMaterial & getMaterialByElement(const Element & element) const; /** * \brief Returns the ith material of *this. * \param i The ith material */ const StructuralMaterial & getMaterial(UInt material_index) const; const StructuralMaterial & getMaterial(const ID & name) const; /** * \brief Returns the number of the different materials inside *this. */ UInt getNbMaterials() const { return materials.size(); } /* ------------------------------------------------------------------------ */ /* Boundaries (structural_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// Compute Linear load function set in global axis template - void computeForcesByGlobalTractionArray(const Array & tractions); + void computeForcesByGlobalTractionArray(const Array & traction_global); /// Compute Linear load function set in local axis template void computeForcesByLocalTractionArray(const Array & tractions); /// compute force vector from a function(x,y,momentum) that describe stresses // template // void computeForcesFromFunction(BoundaryFunction in_function, // BoundaryFunctionType function_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array std::unique_ptr> displacement_rotation; /// velocities array std::unique_ptr> velocity; /// accelerations array std::unique_ptr> acceleration; /// forces array std::unique_ptr> internal_force; /// forces array std::unique_ptr> external_force; /// lumped mass array std::unique_ptr> mass; /// boundaries array std::unique_ptr> blocked_dofs; /// stress array ElementTypeMapArray stress; ElementTypeMapArray element_material; // Define sets of beams ElementTypeMapArray set_ID; /// number of degre of freedom UInt nb_degree_of_freedom; // Rotation matrix ElementTypeMapArray rotation_matrix; // /// analysis method check the list in akantu::AnalysisMethod // AnalysisMethod method; /// flag defining if the increment must be computed or not bool increment_flag; bool need_to_reassemble_mass{true}; bool need_to_reassemble_stiffness{true}; /* ------------------------------------------------------------------------ */ std::vector materials; std::map materials_names_to_id; }; } // namespace akantu #include "structural_mechanics_model_inline_impl.hh" #endif /* AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ */