diff --git a/examples/python/solver_callback/solver_callback.py b/examples/python/solver_callback/solver_callback.py index 3517ece2a..85729bc31 100644 --- a/examples/python/solver_callback/solver_callback.py +++ b/examples/python/solver_callback/solver_callback.py @@ -1,172 +1,177 @@ #!/usr/bin/env python3 """ solver_callback.py: solver_callback overload example""" __author__ = "Nicolas Richart" __credits__ = [ "Guillaume Anciaux ", "Nicolas Richart ", ] __copyright__ = "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale" \ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ " en Mécanique des Solides)" __license__ = "LGPLv3" import numpy as np import akantu as aka class SolverCallback(aka.InterceptSolverCallback): def __init__(self, model): - super().__init__(model.getDOFManager()) + super().__init__(model) self.model = model mesh = model.getMesh() left = mesh.getElementGroup("Left").getNodeGroup().getNodes() right = mesh.getElementGroup("Right").getNodeGroup().getNodes() position = mesh.getNodes() self.pair = [] for node_l in left: node_l = int(node_l) for node_r in right: node_r = int(node_r) if abs(position[node_r, 1] - position[node_l, 1]) < 1e-6: self.pair.append([node_l, node_r]) blocked_dofs = model.getBlockedDOFs() self.periodic_K_modif = aka.TermsToAssemble("displacement", "displacement") matrix_type = self.model.getMatrixType("K") for p in self.pair: #blocked_dofs[p[1]] = True + # a u_{i, x} + b u_{j, x} = 0 + # self.periodic_K_modif(i*dim + aka._x, i*dim + aka._x, a) + # self.periodic_K_modif(i*dim + aka._x, j*dim + aka._x, b) + self.periodic_K_modif(p[0]*2, p[0]*2, 1) self.periodic_K_modif(p[0]*2, p[1]*2, -1) if matrix_type == aka._unsymmetric: self.periodic_K_modif(p[1]*2, p[0]*2, -1) self.periodic_K_modif(p[1]*2, p[1]*2, 1) self.first = True self.k_release = -1 def assembleMatrix(self, matrix_id): self.model.assembleMatrix(matrix_id) if matrix_id == "K": release = self.model.getDOFManager().getMatrix("K").getRelease() if release == self.k_release: return if self.first: self.model.getDOFManager().getMatrix("K").saveMatrix("K0.mtx") self.model.getDOFManager().assemblePreassembledMatrix( "K", self.periodic_K_modif) + if self.first: self.model.getDOFManager().getMatrix("K").saveMatrix("K1.mtx") self.k_release = self.model.getDOFManager().getMatrix("K").getRelease() self.first = False def assembleResidual(self): displacement = self.model.getDisplacement() force = np.zeros(displacement.shape) for p in self.pair: force[p[0], 0] += displacement[p[0], 0] - displacement[p[1], 0] force[p[1], 0] += displacement[p[1], 0] - displacement[p[0], 0] self.model.getDOFManager().assembleToResidual('displacement', force, -1.); self.model.assembleResidual(); # ----------------------------------------------------------------------------- def main(): spatial_dimension = 2 mesh_file = 'bar.msh' max_steps = 250 time_step = 1e-3 aka.parseInput('material.dat') # ------------------------------------------------------------------------- # Initialization # ------------------------------------------------------------------------- mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) model = aka.SolidMechanicsModel(mesh) model.initFull(_analysis_method=aka._implicit_dynamic) model.setBaseName("solver_callback") model.addDumpFieldVector("displacement") model.addDumpFieldVector("acceleration") model.addDumpFieldVector("velocity") model.addDumpFieldVector("internal_force") model.addDumpFieldVector("external_force") model.addDumpField("strain") model.addDumpField("stress") model.addDumpField("blocked_dofs") # ------------------------------------------------------------------------- # boundary conditions # ------------------------------------------------------------------------- model.applyBC(aka.FixedValue(0, aka._y), "YBlocked") # ------------------------------------------------------------------------- # initial conditions # ------------------------------------------------------------------------- displacement = model.getDisplacement() velocity = model.getVelocity() nb_nodes = mesh.getNbNodes() position = mesh.getNodes() L = 1 # pulse_width A = 0.01 v = np.sqrt(model.getMaterial(0).getReal('E') / model.getMaterial(0).getReal('rho')) k = 0.1 * 2 * np.pi * 3 / L t = 0. velocity[:, 0] = k * v * A * np.sin(k * ((position[:, 0] - 5.) - v * t)) displacement[:, 0] = A * np.cos(k * ((position[:, 0] - 5.) - v * t)) # ------------------------------------------------------------------------- # timestep value computation # ------------------------------------------------------------------------- time_factor = 0.8 stable_time_step = model.getStableTimeStep() * time_factor print("Stable Time Step = {0}".format(stable_time_step)) print("Required Time Step = {0}".format(time_step)) time_step = stable_time_step * time_factor model.setTimeStep(time_step) solver_callback = SolverCallback(model) solver = model.getNonLinearSolver() solver.set("max_iterations", 100) solver.set("threshold", 1e-7) # ------------------------------------------------------------------------- # loop for evolution of motion dynamics # ------------------------------------------------------------------------- print("step,step * time_step,epot,ekin,epot + ekin") for step in range(0, max_steps + 1): model.solveStep(solver_callback) #model.solveStep() if step % 10 == 0: model.dump() epot = model.getEnergy('potential') ekin = model.getEnergy('kinetic') # output energy calculation to screen print("{0},{1},{2},{3},{4}".format(step, step * time_step, epot, ekin, (epot + ekin))) return # ----------------------------------------------------------------------------- if __name__ == "__main__": main() diff --git a/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index f282187b1..2d9064dc5 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,919 +1,919 @@ /** * @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: Fri Apr 09 2021 * * @brief Implementation of HeatTransferModel class * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_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, std::shared_ptr dof_manager) : Model(mesh, ModelType::_heat_transfer_model, dof_manager, 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->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) { +MatrixType HeatTransferModel::getMatrixType(const ID & matrix_id) const { 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( const std::string & /* field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/, ElementKind /*element_kind*/) { return nullptr; } /* -------------------------------------------------------------------------- */ 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(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ 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 925f26428..af4f77627 100644 --- a/src/model/heat_transfer/heat_transfer_model.hh +++ b/src/model/heat_transfer/heat_transfer_model.hh @@ -1,316 +1,316 @@ /** * @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 Mar 15 2021 * * @brief Model of Heat Transfer * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #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", std::shared_ptr dof_manager = nullptr); ~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; + MatrixType getMatrixType(const ID & matrix_id) const 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; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Density, density, Real); AKANTU_GET_MACRO(Capacity, capacity, Real); /// 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_couplers/coupler_solid_phasefield.cc b/src/model/model_couplers/coupler_solid_phasefield.cc index 3438f6798..2f3810470 100644 --- a/src/model/model_couplers/coupler_solid_phasefield.cc +++ b/src/model/model_couplers/coupler_solid_phasefield.cc @@ -1,630 +1,630 @@ /** * @file coupler_solid_phasefield.cc * * @author Mohit Pundir * * @date creation: Mon Jun 24 2019 * @date last modification: Fri Apr 02 2021 * * @brief class for coupling of solid mechancis and phase model * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "coupler_solid_phasefield.hh" #include "dumpable_inline_impl.hh" #include "element_synchronizer.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidPhaseField::CouplerSolidPhaseField(Mesh & mesh, UInt dim, const ID & id, const ModelType model_type) : Model(mesh, model_type, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("CouplerSolidPhaseField", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_phasefield", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_phasefield", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModel(mesh, Model::spatial_dimension, "solid_mechanics_model"); phase = new PhaseFieldModel(mesh, Model::spatial_dimension, "phase_field_model"); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_csp_damage); this->registerSynchronizer(synchronizer, SynchronizationTag::_csp_strain); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidPhaseField::~CouplerSolidPhaseField() = default; /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); solid->initFull(_analysis_method = this->method); phase->initFull(_analysis_method = this->method); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidPhaseField::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); auto & phase_model_solver = aka::as_type(*phase); phase_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); } /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidPhaseField::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } 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); } } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidPhaseField::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidPhaseField::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_linear; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleResidual() { // computes the internal forces this->assembleInternalForces(); auto & solid_internal_force = solid->getInternalForce(); auto & solid_external_force = solid->getExternalForce(); auto & phasefield_internal_force = phase->getInternalForce(); auto & phasefield_external_force = phase->getExternalForce(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", solid_external_force, 1); this->getDOFManager().assembleToResidual("displacement", solid_internal_force, 1); this->getDOFManager().assembleToResidual("damage", phasefield_external_force, 1); this->getDOFManager().assembleToResidual("damage", phasefield_internal_force, 1); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); auto & solid_internal_force = solid->getInternalForce(); auto & solid_external_force = solid->getExternalForce(); auto & phasefield_internal_force = phase->getInternalForce(); auto & phasefield_external_force = phase->getExternalForce(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", solid_external_force, 1); this->getDOFManager().assembleToResidual("displacement", solid_internal_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("damage", phasefield_external_force, 1); this->getDOFManager().assembleToResidual("damage", phasefield_internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::predictor() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.predictor(); auto & phase_model_solver = aka::as_type(*phase); phase_model_solver.predictor(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::corrector() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.corrector(); auto & phase_model_solver = aka::as_type(*phase); phase_model_solver.corrector(); } /* -------------------------------------------------------------------------- */ -MatrixType CouplerSolidPhaseField::getMatrixType(const ID & matrix_id) { +MatrixType CouplerSolidPhaseField::getMatrixType(const ID & matrix_id) const { if (matrix_id == "K") { return _symmetric; } if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::beforeSolveStep() { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.beforeSolveStep(); auto & phase_solver_callback = aka::as_type(*phase); phase_solver_callback.beforeSolveStep(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::afterSolveStep(bool converged) { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.afterSolveStep(converged); auto & phase_solver_callback = aka::as_type(*phase); phase_solver_callback.afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); phase->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); phase->assembleStiffnessMatrix(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMassLumped(GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMass(GhostType ghost_type) { solid->assembleMass(ghost_type); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::computeDamageOnQuadPoints( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = phase->getFEEngine(); auto & mesh = phase->getMesh(); auto nb_materials = solid->getNbMaterials(); auto nb_phasefields = phase->getNbPhaseFields(); AKANTU_DEBUG_ASSERT( nb_phasefields == nb_materials, "The number of phasefields and materials should be equal"); for (auto index : arange(nb_materials)) { auto & material = solid->getMaterial(index); for (auto index2 : arange(nb_phasefields)) { auto & phasefield = phase->getPhaseField(index2); if (phasefield.getName() == material.getName()) { switch (spatial_dimension) { case 1: { auto & mat = static_cast &>(material); auto & damage = mat.getDamage(); for (const auto & type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase->getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } case 2: { auto & mat = static_cast &>(material); auto & damage = mat.getDamage(); for (const auto & type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase->getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } default: auto & mat = static_cast &>(material); auto & damage = mat.getDamage(); for (const auto & type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase->getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } } } } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::computeStrainOnQuadPoints( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & mesh = solid->getMesh(); auto nb_materials = solid->getNbMaterials(); auto nb_phasefields = phase->getNbPhaseFields(); AKANTU_DEBUG_ASSERT( nb_phasefields == nb_materials, "The number of phasefields and materials should be equal"); for (auto index : arange(nb_materials)) { auto & material = solid->getMaterial(index); for (auto index2 : arange(nb_phasefields)) { auto & phasefield = phase->getPhaseField(index2); if (phasefield.getName() == material.getName()) { auto & strain_on_qpoints = phasefield.getStrain(); const auto & gradu_on_qpoints = material.getGradU(); for (const auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { auto & strain_on_qpoints_vect = strain_on_qpoints(type, ghost_type); const auto & gradu_on_qpoints_vect = gradu_on_qpoints(type, ghost_type); for (auto && values : zip(make_view(strain_on_qpoints_vect, spatial_dimension, spatial_dimension), make_view(gradu_on_qpoints_vect, spatial_dimension, spatial_dimension))) { auto & strain = std::get<0>(values); const auto & grad_u = std::get<1>(values); gradUToEpsilon(grad_u, strain); } } break; } } } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::solve(const ID & solid_solver_id, const ID & phase_solver_id) { solid->solveStep(solid_solver_id); this->computeStrainOnQuadPoints(_not_ghost); phase->solveStep(phase_solver_id); this->computeDamageOnQuadPoints(_not_ghost); solid->assembleInternalForces(); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { for (UInt j = 0; j < Model::spatial_dimension; ++j) { epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i)); } } } /* ------------------------------------------------------------------------- */ bool CouplerSolidPhaseField::checkConvergence(Array & u_new, Array & u_old, Array & d_new, Array & d_old) { const Array & blocked_dofs = solid->getBlockedDOFs(); UInt nb_degree_of_freedom = u_new.size(); auto u_n_it = u_new.begin(); auto u_o_it = u_old.begin(); auto bld_it = blocked_dofs.begin(); Real norm = 0; for (UInt n = 0; n < nb_degree_of_freedom; ++n, ++u_n_it, ++u_o_it, ++bld_it) { if ((!*bld_it)) { norm += (*u_n_it - *u_o_it) * (*u_n_it - *u_o_it); } } norm = std::sqrt(norm); auto d_n_it = d_new.begin(); auto d_o_it = d_old.begin(); nb_degree_of_freedom = d_new.size(); Real norm2 = 0; for (UInt i = 0; i < nb_degree_of_freedom; ++i) { norm2 += (*d_n_it - *d_o_it); } norm2 = std::sqrt(norm2); Real error = std::max(norm, norm2); Real tolerance = 1e-8; return error < tolerance; } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) { return solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { return solid->createNodalFieldReal(field_name, group_name, padding_flag); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { return solid->createNodalFieldBool(field_name, group_name, padding_flag); std::shared_ptr field; return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createElementalField( const std::string &, const std::string &, bool, UInt, ElementKind) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidPhaseField::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -----------------------------------------------------------------------*/ void CouplerSolidPhaseField::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* ------------------------------------------------------------------------*/ void CouplerSolidPhaseField::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ----------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_phasefield.hh b/src/model/model_couplers/coupler_solid_phasefield.hh index 38503760a..fa64644dd 100644 --- a/src/model/model_couplers/coupler_solid_phasefield.hh +++ b/src/model/model_couplers/coupler_solid_phasefield.hh @@ -1,294 +1,294 @@ /** * @file coupler_solid_phasefield.hh * * @author Mohit Pundir * * @date creation: Mon Jun 24 2019 * @date last modification: Wed Jun 23 2021 * * @brief class for coupling of solid mechancis and phasefield model * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "material.hh" #include "material_phasefield.hh" #include "model.hh" #include "phase_field_model.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ #define __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ /* ------------------------------------------------------------------------ */ /* Coupling : Solid Mechanics / PhaseField */ /* ------------------------------------------------------------------------ */ namespace akantu { template class IntegratorGauss; template class ShapeLagrange; class DOFManager; } // namespace akantu namespace akantu { class CouplerSolidPhaseField : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: CouplerSolidPhaseField( Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "coupler_solid_phasefield", ModelType model_type = ModelType::_coupler_solid_phasefield); ~CouplerSolidPhaseField() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize the complete model void initFullImpl(const ModelOptions & options) override; /// initialize the modelType void initModel() override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); public: /// computes damage on quad points for solid mechanics model from /// damage array from phasefield model void computeDamageOnQuadPoints(const GhostType & /*ghost_type*/); /// computes strain on quadrature points for phasefield model from /// displacement gradient from solid mechanics model void computeStrainOnQuadPoints(const GhostType & ghost_type); /// solve the coupled model void solve(const ID & solid_solver_id = "", const ID & phase_solver_id = ""); private: /// computes small strain from displacement gradient void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); /// test the convergence criteria bool checkConvergence(Array & /*u_new*/, Array & /*u_old*/, Array & /*d_new*/, Array & /*d_old*/); 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; } + bool canSplitResidual() const override { return true; } /// get the type of matrix needed - MatrixType getMatrixType(const ID & matrix_id) override; + MatrixType getMatrixType(const ID & matrix_id) const 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 model to instantiate the matricess when needed void initSolver(TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType /*non_linear_solver_type*/) 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; /// solve the coupled model // void solveStep(const ID & solver_id = "") override; /// solve a step using a given pre instantiated time step solver and /// non linear solver with a user defined callback instead of the /// model itself /!\ This can mess up everything // void solveStep(SolverCallback & callback, const ID & solver_id = "") // override; /* ------------------------------------------------------------------------ */ /* Mass matrix for solid mechanics model */ /* ------------------------------------------------------------------------ */ 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); protected: /* -------------------------------------------------------------------------- */ TimeStepSolverType getDefaultSolverType() const override; /* -------------------------------------------------------------------------- */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass; } /* ------------------------------------------------------------------------ */ public: // DataAccessor UInt getNbData(const Array & /*elements*/, const SynchronizationTag & /*tag*/) const override { return 0; } void packData(CommunicationBuffer & /*buffer*/, const Array & /*element*/, const SynchronizationTag & /*tag*/) const override {} void unpackData(CommunicationBuffer & /*buffer*/, const Array & /*element*/, const SynchronizationTag & /*tag*/) override {} UInt getNbData(__attribute__((unused)) const Array & indexes, __attribute__((unused)) const SynchronizationTag & tag) const override { return 0; } void packData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Array & dofs, __attribute__((unused)) const SynchronizationTag & tag) const override {} void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Array & dofs, __attribute__((unused)) const SynchronizationTag & tag) override {} /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the solid mechanics model AKANTU_GET_MACRO(SolidMechanicsModel, *solid, SolidMechanicsModel &); /// get the contact mechanics model AKANTU_GET_MACRO(PhaseFieldModel, *phase, PhaseFieldModel &); /* ------------------------------------------------------------------------ */ /* 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; 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; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ private: /// solid mechanics model SolidMechanicsModel * solid{nullptr}; /// phasefield model PhaseFieldModel * phase{nullptr}; Array * displacement{nullptr}; /// Array * displacement_increment{nullptr}; /// external forces array Array * external_force{nullptr}; }; } // namespace akantu #endif /* __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ */ diff --git a/src/model/phase_field/phase_field_model.cc b/src/model/phase_field/phase_field_model.cc index fadfd6ffe..84e0a8418 100644 --- a/src/model/phase_field/phase_field_model.cc +++ b/src/model/phase_field/phase_field_model.cc @@ -1,635 +1,635 @@ /** * @file phase_field_model.cc * * @author Mohit Pundir * * @date creation: Tue Sep 04 2018 * @date last modification: Wed Jun 23 2021 * * @brief Implementation of PhaseFieldModel class * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "phase_field_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 { /* -------------------------------------------------------------------------- */ PhaseFieldModel::PhaseFieldModel(Mesh & mesh, UInt dim, const ID & id, ModelType model_type) : Model(mesh, model_type, dim, id), phasefield_index("phasefield index", id), phasefield_local_numbering("phasefield local numbering", id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("PhaseFieldFEEngine", mesh, Model::spatial_dimension); #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("phase_field", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif // AKANTU_USE_IOHELPER phasefield_selector = std::make_shared(phasefield_index); this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_damage); this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_driving); this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_history); this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_energy); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseFieldModel::~PhaseFieldModel() = default; /* -------------------------------------------------------------------------- */ -MatrixType PhaseFieldModel::getMatrixType(const ID & matrix_id) { +MatrixType PhaseFieldModel::getMatrixType(const ID & matrix_id) const { if (matrix_id == "K") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initModel() { auto & fem = this->getFEEngine(); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initFullImpl(const ModelOptions & options) { phasefield_index.initialize(mesh, _element_kind = _ek_not_defined, _default_value = UInt(-1), _with_nb_element = true); phasefield_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true); Model::initFullImpl(options); // initialize the phasefields if (!this->parser.getLastParsedFile().empty()) { this->instantiatePhaseFields(); this->initPhaseFields(); } this->initBC(*this, *damage, *external_force); } /* -------------------------------------------------------------------------- */ PhaseField & PhaseFieldModel::registerNewPhaseField(const ParserSection & section) { std::string phase_name; std::string phase_type = section.getName(); std::string opt_param = section.getOption(); try { std::string tmp = section.getParameter("name"); phase_name = tmp; /** this can seam weird, but there is an ambiguous * operator overload that i couldn't solve. @todo remove * the weirdness of this code */ } catch (debug::Exception &) { AKANTU_ERROR("A phasefield of type \'" << phase_type << "\' in the input file has been defined without a name!"); } PhaseField & phase = this->registerNewPhaseField(phase_name, phase_type, opt_param); phase.parseSection(section); return phase; } /* -------------------------------------------------------------------------- */ PhaseField & PhaseFieldModel::registerNewPhaseField(const ID & phase_name, const ID & phase_type, const ID & opt_param) { AKANTU_DEBUG_ASSERT(phasefields_names_to_id.find(phase_name) == phasefields_names_to_id.end(), "A phasefield with this name '" << phase_name << "' has already been registered. " << "Please use unique names for phasefields"); UInt phase_count = phasefields.size(); phasefields_names_to_id[phase_name] = phase_count; std::stringstream sstr_phase; sstr_phase << this->id << ":" << phase_count << ":" << phase_type; ID mat_id = sstr_phase.str(); std::unique_ptr phase = PhaseFieldFactory::getInstance().allocate( phase_type, opt_param, *this, mat_id); phasefields.push_back(std::move(phase)); return *(phasefields.back()); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::instantiatePhaseFields() { ParserSection model_section; bool is_empty; std::tie(model_section, is_empty) = this->getParserSection(); if (not is_empty) { auto model_phasefields = model_section.getSubSections(ParserType::_phasefield); for (const auto & section : model_phasefields) { this->registerNewPhaseField(section); } } auto sub_sections = this->parser.getSubSections(ParserType::_phasefield); for (const auto & section : sub_sections) { this->registerNewPhaseField(section); } if (phasefields.empty()) { AKANTU_EXCEPTION("No phasefields where instantiated for the model" << getID()); } are_phasefields_instantiated = true; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initPhaseFields() { AKANTU_DEBUG_ASSERT(phasefields.size() != 0, "No phasefield to initialize !"); if (!are_phasefields_instantiated) { instantiatePhaseFields(); } this->assignPhaseFieldToElements(); for (auto & phasefield : phasefields) { /// init internals properties phasefield->initPhaseField(); } this->synchronize(SynchronizationTag::_smm_init_mat); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assignPhaseFieldToElements( const ElementTypeMapArray * filter) { for_each_element( mesh, [&](auto && element) { UInt phase_index = (*phasefield_selector)(element); AKANTU_DEBUG_ASSERT( phase_index < phasefields.size(), "The phasefield selector returned an index that does not exists"); phasefield_index(element) = phase_index; }, _element_filter = filter, _ghost_type = _not_ghost); for_each_element( mesh, [&](auto && element) { auto phase_index = phasefield_index(element); auto index = phasefields[phase_index]->addElement(element); phasefield_local_numbering(element) = index; }, _element_filter = filter, _ghost_type = _not_ghost); // synchronize the element phasefield arrays this->synchronize(SynchronizationTag::_material_id); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else { AKANTU_ERROR("Unknown Matrix ID for PhaseFieldModel : " << matrix_id); } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::predictor() { // AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::corrector() { // AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->damage, 1, "damage"); this->allocNodalField(this->external_force, 1, "external_force"); this->allocNodalField(this->internal_force, 1, "internal_force"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); this->allocNodalField(this->previous_damage, 1, "previous_damage"); if (!dof_manager.hasDOFs("damage")) { dof_manager.registerDOFs("damage", *this->damage, _dst_nodal); dof_manager.registerBlockedDOFs("damage", *this->blocked_dofs); dof_manager.registerDOFsPrevious("damage", *this->previous_damage); } if (time_step_solver_type == TimeStepSolverType::_dynamic) { AKANTU_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ FEEngine & PhaseFieldModel::getFEEngineBoundary(const ID & name) { return dynamic_cast(getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ std::tuple PhaseFieldModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } 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 PhaseFieldModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["damage"] = IntegrationSchemeType::_central_difference; options.solution_type["damage"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_linear; options.integration_scheme_type["damage"] = IntegrationSchemeType::_pseudo_time; options.solution_type["damage"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["damage"] = IntegrationSchemeType::_backward_euler; options.solution_type["damage"] = IntegrationScheme::_damage; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::beforeSolveStep() { for (auto & phasefield : phasefields) { phasefield->beforeSolveStep(); } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::afterSolveStep(bool converged) { if (not converged) { return; } for (auto && values : zip(*damage, *previous_damage)) { auto & dam = std::get<0>(values); auto & prev_dam = std::get<1>(values); dam -= prev_dam; prev_dam = dam; } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleStiffnessMatrix() { AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().zeroMatrix("K"); for (auto & phasefield : phasefields) { phasefield->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleResidual() { this->assembleInternalForces(); this->getDOFManager().assembleToResidual("damage", *this->external_force, 1); this->getDOFManager().assembleToResidual("damage", *this->internal_force, 1); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleInternalForces() { AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->zero(); // communicate the driving forces AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_pfm_driving); // assemble the forces due to local driving forces AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & phasefield : phasefields) { phasefield->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant driving forces"); this->waitEndSynchronize(SynchronizationTag::_pfm_driving); // assemble the residual due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleLumpedMatrix(const ID & /*matrix_id*/) {} /* -------------------------------------------------------------------------- */ void PhaseFieldModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("phase_field").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case SynchronizationTag::_pfm_damage: { size += nb_nodes_per_element * sizeof(Real); // damage break; } case SynchronizationTag::_pfm_driving: { size += getNbIntegrationPoints(elements) * sizeof(Real); break; } case SynchronizationTag::_pfm_history: { size += getNbIntegrationPoints(elements) * sizeof(Real); break; } case SynchronizationTag::_pfm_energy: { size += getNbIntegrationPoints(elements) * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } return size; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::packData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Array & elements, __attribute__((unused)) const SynchronizationTag & tag) const {} /* -------------------------------------------------------------------------- */ void PhaseFieldModel::unpackData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Array & elements, __attribute__((unused)) const SynchronizationTag & tag) {} /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & indexes, const SynchronizationTag & tag) const { UInt size = 0; UInt nb_nodes = indexes.size(); switch (tag) { case SynchronizationTag::_pfm_damage: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } return size; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const { for (auto index : indexes) { switch (tag) { case SynchronizationTag::_pfm_damage: { buffer << (*damage)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) { for (auto index : indexes) { switch (tag) { case SynchronizationTag::_pfm_damage: { buffer >> (*damage)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER std::shared_ptr PhaseFieldModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool /*unused*/) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs.get(); return mesh.createNodalField(uint_nodal_fields[field_name], group_name); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool /*unused*/) { std::map *> real_nodal_fields; real_nodal_fields["damage"] = damage.get(); real_nodal_fields["external_force"] = external_force.get(); real_nodal_fields["internal_force"] = internal_force.get(); return mesh.createNodalField(real_nodal_fields[field_name], group_name); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createElementalField( const std::string & field_name, const std::string & group_name, bool /*unused*/, UInt /*unused*/, ElementKind element_kind) { if (field_name == "partitions") { return mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); } std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createElementalField(const std::string &, const std::string &, bool, const UInt &, ElementKind) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void PhaseFieldModel::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Phase Field Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; damage->printself(stream, indent + 2); external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + phasefield information [" << std::endl; stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } } // namespace akantu diff --git a/src/model/phase_field/phase_field_model.hh b/src/model/phase_field/phase_field_model.hh index c9c90a530..f14f7eeb4 100644 --- a/src/model/phase_field/phase_field_model.hh +++ b/src/model/phase_field/phase_field_model.hh @@ -1,338 +1,338 @@ /** * @file phase_field_model.hh * * @author Mohit Pundir * * @date creation: Tue Sep 04 2018 * @date last modification: Wed Jun 23 2021 * * @brief Model class for Phase Field problem * * * @section LICENSE * * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASE_FIELD_MODEL_HH__ #define __AKANTU_PHASE_FIELD_MODEL_HH__ namespace akantu { class PhaseField; class PhaseFieldSelector; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class PhaseFieldModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using FEEngineType = FEEngineTemplate; PhaseFieldModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "phase_field_model", ModelType model_type = ModelType::_phase_field_model); ~PhaseFieldModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// generic function to initialize everything ready for explicit dynamics void initFullImpl(const ModelOptions & options) override; /// initialize all internal array for phasefields void initPhaseFields(); /// allocate all vectors void initSolver(TimeStepSolverType /*unused*/, NonLinearSolverType /*unused*/) override; /// initialize the model void initModel() override; /// predictor void predictor() override; /// corrector void corrector() override; /// compute the heat flux void assembleResidual() override; /// get the type of matrix needed - MatrixType getMatrixType(const ID & /*unused*/) override; + MatrixType getMatrixType(const ID & /*unused*/) const override; /// callback to assemble a Matrix void assembleMatrix(const ID & /*unused*/) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & /*unused*/) override; std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ - /* Materials (phase_field_model.cc) */ + /* Materials (phase_field_model.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty phasefield of a given type PhaseField & registerNewPhaseField(const ID & phase_name, const ID & phase_type, const ID & opt_param); /// reassigns phasefields depending on the phasefield selector void reassignPhaseField(); protected: /// register a phasefield in the dynamic database PhaseField & registerNewPhaseField(const ParserSection & phase_section); /// read the phasefield files to instantiate all the phasefields void instantiatePhaseFields(); /// set the element_id_by_phasefield and add the elements to the good /// phasefields void assignPhaseFieldToElements( const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Methods for static */ /* ------------------------------------------------------------------------ */ public: /// assembles the phasefield stiffness matrix virtual void assembleStiffnessMatrix(); /// compute the internal forces virtual void assembleInternalForces(); // compute the internal forces void assembleInternalForces(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Methods for dynamic */ /* ------------------------------------------------------------------------ */ public: /// set the stable timestep void setTimeStep(Real time_step, const ID & solver_id = "") override; protected: /// 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; /* ------------------------------------------------------------------------ */ /* 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 & indexes, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the damage array AKANTU_GET_MACRO_DEREF_PTR(Damage, damage); AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Damage, damage); /// get the PhaseFieldModel::internal_force vector (internal forces) AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force); AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force); /// get the PhaseFieldModel::external_force vector (external forces) AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force); AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force); /// get the PhaseFieldModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the PhaseFieldModel::blocked_dofs vector AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs); /// get an iterable on the phasefields inline decltype(auto) getPhaseFields(); /// get an iterable on the phasefields inline decltype(auto) getPhaseFields() const; /// get a particular phasefield (by phasefield index) inline PhaseField & getPhaseField(UInt mat_index); /// get a particular phasefield (by phasefield index) inline const PhaseField & getPhaseField(UInt mat_index) const; /// get a particular phasefield (by phasefield name) inline PhaseField & getPhaseField(const std::string & name); /// get a particular phasefield (by phasefield name) inline const PhaseField & getPhaseField(const std::string & name) const; /// get a particular phasefield id from is name inline UInt getPhaseFieldIndex(const std::string & name) const; /// give the number of phasefields inline UInt getNbPhaseFields() const { return phasefields.size(); } /// give the phasefield internal index from its id Int getInternalIndexFromID(const ID & id) const; AKANTU_GET_MACRO(PhaseFieldByElement, phasefield_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(PhaseFieldLocalNumbering, phasefield_local_numbering, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PhaseFieldByElement, phasefield_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(PhaseFieldByElement, phasefield_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PhaseFieldLocalNumbering, phasefield_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(PhaseFieldLocalNumbering, phasefield_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(PhaseFieldSelector, *phasefield_selector, PhaseFieldSelector &); AKANTU_SET_MACRO(PhaseFieldSelector, phasefield_selector, std::shared_ptr); FEEngine & getFEEngineBoundary(const ID & name = "") 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; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// number of iterations UInt n_iter; /// damage array std::unique_ptr> damage; /// damage array at the previous time step std::unique_ptr> previous_damage; /// boundary vector std::unique_ptr> blocked_dofs; /// external force vector std::unique_ptr> external_force; /// residuals array std::unique_ptr> internal_force; /// Arrays containing the phasefield index for each element ElementTypeMapArray phasefield_index; /// Arrays containing the position in the element filter of the phasefield /// (phasefield's local numbering) ElementTypeMapArray phasefield_local_numbering; /// class defining of to choose a phasefield std::shared_ptr phasefield_selector; /// mapping between phasefield name and phasefield internal id std::map phasefields_names_to_id; /// list of used phasefields std::vector> phasefields; /// tells if the phasefield are instantiated bool are_phasefields_instantiated{false}; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "parser.hh" #include "phasefield.hh" #include "phase_field_model_inline_impl.cc" /* -------------------------------------------------------------------------- */ #endif diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index 69da5e251..6419ff3b7 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,627 +1,626 @@ /** * @file structural_mechanics_model.cc * * @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: Mon Mar 15 2021 * * @brief Model implementation for Structural Mechanics elements * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model.hh" #include "dof_manager.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "shape_structural.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumpable_inline_impl.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" #include "group_manager_inline_impl.hh" #endif /* -------------------------------------------------------------------------- */ #include "structural_element_bernoulli_beam_2.hh" #include "structural_element_bernoulli_beam_3.hh" #include "structural_element_kirchhoff_shell.hh" /* -------------------------------------------------------------------------- */ //#include "structural_mechanics_model_inline_impl.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt StructuralMechanicsModel::getNbDegreeOfFreedom(ElementType type) { UInt ndof = 0; #define GET_(type) ndof = ElementClass::getNbDegreeOfFreedom() AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_, _ek_structural); #undef GET_ return ndof; } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id) : Model(mesh, ModelType::_structural_mechanics_model, dim, id), f_m2a(1.0), stress("stress", id), element_material("element_material", id), set_ID("beam sets", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("StructuralMechanicsFEEngine", mesh, spatial_dimension); if (spatial_dimension == 2) { nb_degree_of_freedom = 3; } else if (spatial_dimension == 3) { nb_degree_of_freedom = 6; } else { AKANTU_TO_IMPLEMENT(); } #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("structural_mechanics_model", id, true); #endif this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); this->initDOFManager(); this->dumper_default_element_kind = _ek_structural; mesh.getElementalData("extra_normal") .initialize(mesh, _element_kind = _ek_structural, _nb_component = spatial_dimension, _with_nb_element = true, _default_value = 0.); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // Initializing stresses ElementTypeMap stress_components; /// TODO this is ugly af, maybe add a function to FEEngine for (auto && type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { UInt nb_components = 0; // Getting number of components for each element type #define GET_(type) nb_components = ElementClass::getNbStressComponents() AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_); #undef GET_ stress_components(nb_components, type); } stress.initialize( getFEEngine(), _spatial_dimension = _all_dimensions, _element_kind = _ek_structural, _nb_component = [&stress_components](ElementType type, GhostType /*unused*/) -> UInt { return stress_components(type); }); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFEEngineBoundary() { /// TODO: this function should not be reimplemented /// we're just avoiding a call to Model::initFEEngineBoundary() } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { AKANTU_DEBUG_IN(); this->allocNodalField(displacement_rotation, nb_degree_of_freedom, "displacement"); this->allocNodalField(external_force, nb_degree_of_freedom, "external_force"); this->allocNodalField(internal_force, nb_degree_of_freedom, "internal_force"); this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs"); auto & dof_manager = this->getDOFManager(); if (not dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *displacement_rotation, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); } if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(velocity, nb_degree_of_freedom, "velocity"); this->allocNodalField(acceleration, nb_degree_of_freedom, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { element_material.initialize(mesh, _element_kind = _ek_structural, _default_value = 0, _with_nb_element = true); getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); if (not need_to_reassemble_stiffness) { return; } if (not getDOFManager().hasMatrix("K")) { getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().zeroMatrix("K"); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define ASSEMBLE_STIFFNESS_MATRIX(type) assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } need_to_reassemble_stiffness = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define COMPUTE_STRESS_ON_QUAD(type) computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::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(); return mesh.createNodalField(uint_nodal_fields[field_name], group_name); } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { UInt n; if (spatial_dimension == 2) { n = 2; } else { n = 3; } UInt padding_size = 0; if (padding_flag) { padding_size = 3; } if (field_name == "displacement") { return mesh.createStridedNodalField(displacement_rotation.get(), group_name, n, 0, padding_size); } if (field_name == "velocity") { return mesh.createStridedNodalField(velocity.get(), group_name, n, 0, padding_size); } if (field_name == "acceleration") { return mesh.createStridedNodalField(acceleration.get(), group_name, n, 0, padding_size); } if (field_name == "rotation") { return mesh.createStridedNodalField(displacement_rotation.get(), group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "force") { return mesh.createStridedNodalField(external_force.get(), group_name, n, 0, padding_size); } if (field_name == "external_force") { return mesh.createStridedNodalField(external_force.get(), group_name, n, 0, padding_size); } if (field_name == "momentum") { return mesh.createStridedNodalField(external_force.get(), group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "internal_force") { return mesh.createStridedNodalField(internal_force.get(), group_name, n, 0, padding_size); } if (field_name == "internal_momentum") { return mesh.createStridedNodalField(internal_force.get(), group_name, nb_degree_of_freedom - n, n, padding_size); } return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool /*unused*/, UInt spatial_dimension, ElementKind kind) { std::shared_ptr field; if (field_name == "element_index_by_material") { field = mesh.createElementalField( field_name, group_name, spatial_dimension, kind); } if (field_name == "stress") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(stress); field = mesh.createElementalField( stress, group_name, this->spatial_dimension, kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ /* Virtual methods from SolverCallback */ /* -------------------------------------------------------------------------- */ /// get the type of matrix needed -MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) { +MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) const { return _symmetric; } /// callback to assemble a Matrix void StructuralMechanicsModel::assembleMatrix(const ID & id) { if (id == "K") { assembleStiffnessMatrix(); } else if (id == "M") { assembleMassMatrix(); } } /// callback to assemble a lumped Matrix void StructuralMechanicsModel::assembleLumpedMatrix(const ID & /*id*/) {} /// callback to assemble the residual StructuralMechanicsModel::(rhs) void StructuralMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); auto & dof_manager = getDOFManager(); assembleInternalForce(); dof_manager.assembleToResidual("displacement", *external_force, 1); dof_manager.assembleToResidual("displacement", *internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->assembleInternalForce(); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Virtual methods from Model */ /* -------------------------------------------------------------------------- */ /// get some default values for derived classes std::tuple StructuralMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { 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 StructuralMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce() { internal_force->zero(); computeStresses(); for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { assembleInternalForce(type, _not_ghost); // assembleInternalForce(type, _ghost); } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce(ElementType type, GhostType gt) { auto & fem = getFEEngine(); auto & sigma = stress(type, gt); auto ndof = getNbDegreeOfFreedom(type); auto nb_nodes = mesh.getNbNodesPerElement(type); auto ndof_per_elem = ndof * nb_nodes; Array BtSigma(fem.getNbIntegrationPoints(type) * mesh.getNbElement(type), ndof_per_elem, "BtSigma"); fem.computeBtD(sigma, BtSigma, type, gt); Array intBtSigma(0, ndof_per_elem, "intBtSigma"); fem.integrate(BtSigma, intBtSigma, ndof_per_elem, type, gt); getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force, type, gt, -1.); } /* -------------------------------------------------------------------------- */ Real StructuralMechanicsModel::getKineticEnergy() { if (not this->getDOFManager().hasMatrix("M")) { return 0.; } Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Array Mv(nb_nodes, nb_degree_of_freedom); this->getDOFManager().assembleMatMulVectToArray("displacement", "M", *this->velocity, Mv); for (auto && data : zip(arange(nb_nodes), make_view(Mv, nb_degree_of_freedom), make_view(*this->velocity, nb_degree_of_freedom))) { ekin += std::get<2>(data).dot(std::get<1>(data)) * static_cast(mesh.isLocalOrMasterNode(std::get<0>(data))); } mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum); return ekin / 2.; } /* -------------------------------------------------------------------------- */ Real StructuralMechanicsModel::getPotentialEnergy() { Real epot = 0.; UInt nb_nodes = mesh.getNbNodes(); Array Ku(nb_nodes, nb_degree_of_freedom); this->getDOFManager().assembleMatMulVectToArray( "displacement", "K", *this->displacement_rotation, Ku); for (auto && data : zip(arange(nb_nodes), make_view(Ku, nb_degree_of_freedom), make_view(*this->displacement_rotation, nb_degree_of_freedom))) { epot += std::get<2>(data).dot(std::get<1>(data)) * static_cast(mesh.isLocalOrMasterNode(std::get<0>(data))); } mesh.getCommunicator().allReduce(epot, SynchronizerOperation::_sum); return epot / 2.; } /* -------------------------------------------------------------------------- */ Real StructuralMechanicsModel::getEnergy(const ID & energy) { if (energy == "kinetic") { return getKineticEnergy(); } if (energy == "potential") { return getPotentialEnergy(); } return 0; } -/* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeForcesByLocalTractionArray( const Array & tractions, ElementType type) { AKANTU_DEBUG_IN(); auto nb_element = getFEEngine().getMesh().getNbElement(type); auto nb_nodes_per_element = getFEEngine().getMesh().getNbNodesPerElement(type); auto nb_quad = getFEEngine().getNbIntegrationPoints(type); // check dimension match AKANTU_DEBUG_ASSERT( Mesh::getSpatialDimension(type) == getFEEngine().getElementDimension(), "element type dimension does not match the dimension of boundaries : " << getFEEngine().getElementDimension() << " != " << Mesh::getSpatialDimension(type)); // check size of the vector AKANTU_DEBUG_ASSERT( tractions.size() == nb_quad * nb_element, "the size of the vector should be the total number of quadrature points"); // check number of components AKANTU_DEBUG_ASSERT(tractions.getNbComponent() == nb_degree_of_freedom, "the number of components should be the spatial " "dimension of the problem"); Array Ntbs(nb_element * nb_quad, nb_degree_of_freedom * nb_nodes_per_element); auto & fem = getFEEngine(); fem.computeNtb(tractions, Ntbs, type); // allocate the vector that will contain the integrated values auto name = id + std::to_string(type) + ":integral_boundary"; Array int_funct(nb_element, nb_degree_of_freedom * nb_nodes_per_element, name); // do the integration getFEEngine().integrate(Ntbs, int_funct, nb_degree_of_freedom * nb_nodes_per_element, type); // assemble the result into force vector getDOFManager().assembleElementalArrayLocalArray(int_funct, *external_force, type, _not_ghost, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeForcesByGlobalTractionArray( const Array & traction_global, ElementType type) { AKANTU_DEBUG_IN(); UInt nb_element = mesh.getNbElement(type); UInt nb_quad = getFEEngine().getNbIntegrationPoints(type); Array traction_local(nb_element * nb_quad, nb_degree_of_freedom, id + ":structuralmechanics:imposed_linear_load"); auto R_it = getFEEngineClass() .getShapeFunctions() .getRotations(type) .begin(nb_degree_of_freedom, nb_degree_of_freedom); auto Te_it = traction_global.begin(nb_degree_of_freedom); auto te_it = traction_local.begin(nb_degree_of_freedom); for (UInt e = 0; e < nb_element; ++e, ++R_it) { for (UInt q = 0; q < nb_quad; ++q, ++Te_it, ++te_it) { // turn the traction in the local referential te_it->template mul(*R_it, *Te_it); } } computeForcesByLocalTractionArray(traction_local, type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::afterSolveStep(bool converged) { if (converged) { assembleInternalForce(); } } } // namespace akantu diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh index 770ebb5a0..b3ac31d3f 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.hh +++ b/src/model/structural_mechanics/structural_mechanics_model.hh @@ -1,335 +1,335 @@ /** * @file structural_mechanics_model.hh * * @author Fabian Barras * @author Lucas Frerot * @author Sébastien Hartmann * @author Philip Mueller * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Thu Apr 01 2021 * * @brief Particular implementation of the structural elements in the * StructuralMechanicsModel * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "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; + MatrixType getMatrixType(const ID & matrix_id) const 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 true; } + bool canSplitResidual() const override { return true; } void afterSolveStep(bool converged) override; /// 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; /// 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 void computeForcesByGlobalTractionArray(const Array & traction_global, ElementType type); /// Compute Linear load function set in local axis void computeForcesByLocalTractionArray(const Array & tractions, ElementType 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_ */ diff --git a/test/test_model/test_common/test_model_solver/test_dof_manager_default.cc b/test/test_model/test_common/test_model_solver/test_dof_manager_default.cc index c814910f3..a4ef18cc1 100644 --- a/test/test_model/test_common/test_model_solver/test_dof_manager_default.cc +++ b/test/test_model/test_common/test_model_solver/test_dof_manager_default.cc @@ -1,131 +1,133 @@ /** * @file test_dof_manager_default.cc * * @author Nicolas Richart * * @date creation: Fri Feb 26 2016 * @date last modification: Wed Jan 30 2019 * * @brief Test default dof manager * * * @section LICENSE * * Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager_default.hh" #include "solver_callback.hh" #include "sparse_matrix_aij.hh" #include "time_step_solver.hh" using namespace akantu; /** * =\o-----o-----o-> F * | | * |---- L ----| */ class MySolverCallback : public SolverCallback { public: MySolverCallback(Real F, DOFManagerDefault & dof_manager, UInt nb_dofs = 3) : dof_manager(dof_manager), dispacement(nb_dofs, 1, "disp"), blocked(nb_dofs, 1), forces(nb_dofs, 1), nb_dofs(nb_dofs) { dof_manager.registerDOFs("disp", dispacement, _dst_generic); dof_manager.registerBlockedDOFs("disp", blocked); dispacement.set(0.); forces.set(0.); blocked.set(false); forces(nb_dofs - 1, _x) = F; blocked(0, _x) = true; } - void assembleMatrix(const ID & matrix_id) { + void assembleMatrix(const ID & matrix_id) override { if (matrix_id != "K") return; auto & K = dynamic_cast(dof_manager.getMatrix("K")); K.zero(); for (UInt i = 1; i < nb_dofs - 1; ++i) K.add(i, i, 2.); for (UInt i = 0; i < nb_dofs - 1; ++i) K.add(i, i + 1, -1.); K.add(0, 0, 1); K.add(nb_dofs - 1, nb_dofs - 1, 1); // K *= 1 / L_{el} K *= nb_dofs - 1; } - MatrixType getMatrixType(const ID & matrix_id) { + MatrixType getMatrixType(const ID & matrix_id) const override { if (matrix_id == "K") return _symmetric; return _mt_not_defined; } - void assembleLumpedMatrix(const ID &) {} + void assembleLumpedMatrix(const ID &) override {} - void assembleResidual() { dof_manager.assembleToResidual("disp", forces); } + void assembleResidual() override { + dof_manager.assembleToResidual("disp", forces); + } - void predictor() {} - void corrector() {} + void predictor() override {} + void corrector() override {} DOFManagerDefault & dof_manager; Array dispacement; Array blocked; Array forces; UInt nb_dofs; }; int main(int argc, char * argv[]) { initialize(argc, argv); DOFManagerDefault dof_manager("test_dof_manager"); MySolverCallback callback(10., dof_manager, 11); NonLinearSolver & nls = dof_manager.getNewNonLinearSolver("my_nls", NonLinearSolverType::_linear); TimeStepSolver & tss = dof_manager.getNewTimeStepSolver( "my_tss", TimeStepSolverType::_static, nls, callback); tss.setIntegrationScheme("disp", IntegrationSchemeType::_pseudo_time); tss.solveStep(callback); dof_manager.getMatrix("K").saveMatrix("K_dof_manager_default.mtx"); Array::const_scalar_iterator disp_it = callback.dispacement.begin(); Array::const_scalar_iterator force_it = callback.forces.begin(); Array::const_scalar_iterator blocked_it = callback.blocked.begin(); std::cout << std::setw(8) << "disp" << " " << std::setw(8) << "force" << " " << std::setw(8) << "blocked" << std::endl; for (; disp_it != callback.dispacement.end(); ++disp_it, ++force_it, ++blocked_it) { std::cout << std::setw(8) << *disp_it << " " << std::setw(8) << *force_it << " " << std::setw(8) << std::boolalpha << *blocked_it << std::endl; } finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh b/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh index b1f750ae0..09848ab42 100644 --- a/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh +++ b/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh @@ -1,453 +1,453 @@ /** * @file test_model_solver_my_model.hh * * @author Nicolas Richart * * @date creation: Wed Apr 13 2016 * @date last modification: Fri Jun 26 2020 * * @brief Test default dof manager * * * @section LICENSE * * Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" #include "boundary_condition.hh" #include "communicator.hh" #include "data_accessor.hh" #include "dof_manager_default.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "model_solver.hh" #include "periodic_node_synchronizer.hh" #include "solver_vector_default.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { #ifndef AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH_ #define AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH_ /** * =\o-----o-----o-> F * | | * |---- L ----| */ class MyModel : public ModelSolver, public BoundaryCondition, public DataAccessor { public: MyModel(Real F, Mesh & mesh, bool lumped, const ID & dof_manager_type = "default") : ModelSolver(mesh, ModelType::_model, "model_solver"), nb_dofs(mesh.getNbNodes()), nb_elements(mesh.getNbElement(_segment_2)), lumped(lumped), E(1.), A(1.), rho(1.), mesh(mesh), displacement(nb_dofs, 1, "disp"), velocity(nb_dofs, 1, "velo"), acceleration(nb_dofs, 1, "accel"), blocked(nb_dofs, 1, "blocked"), forces(nb_dofs, 1, "force_ext"), internal_forces(nb_dofs, 1, "force_int"), stresses(nb_elements, 1, "stress"), strains(nb_elements, 1, "strain"), initial_lengths(nb_elements, 1, "L0") { this->initDOFManager(dof_manager_type); this->initBC(*this, displacement, forces); this->getDOFManager().registerDOFs("disp", displacement, _dst_nodal); this->getDOFManager().registerDOFsDerivative("disp", 1, velocity); this->getDOFManager().registerDOFsDerivative("disp", 2, acceleration); this->getDOFManager().registerBlockedDOFs("disp", blocked); displacement.set(0.); velocity.set(0.); acceleration.set(0.); forces.set(0.); blocked.set(false); UInt global_nb_nodes = mesh.getNbGlobalNodes(); for (auto && n : arange(nb_dofs)) { auto global_id = mesh.getNodeGlobalId(n); if (global_id == (global_nb_nodes - 1)) forces(n, _x) = F; if (global_id == 0) blocked(n, _x) = true; } auto cit = this->mesh.getConnectivity(_segment_2).begin(2); auto cend = this->mesh.getConnectivity(_segment_2).end(2); auto L_it = this->initial_lengths.begin(); for (; cit != cend; ++cit, ++L_it) { const Vector & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); *L_it = std::abs(p2 - p1); } this->registerDataAccessor(*this); this->registerSynchronizer( const_cast(this->mesh.getElementSynchronizer()), SynchronizationTag::_user_1); } void assembleLumpedMass() { auto & M = this->getDOFManager().getLumpedMatrix("M"); M.zero(); this->assembleLumpedMass(_not_ghost); if (this->mesh.getNbElement(_segment_2, _ghost) > 0) this->assembleLumpedMass(_ghost); is_lumped_mass_assembled = true; } void assembleLumpedMass(GhostType ghost_type) { Array M(nb_dofs, 1, 0.); Array m_all_el(this->mesh.getNbElement(_segment_2, ghost_type), 2); for (auto && data : zip(make_view(this->mesh.getConnectivity(_segment_2), 2), make_view(m_all_el, 2))) { const auto & conn = std::get<0>(data); auto & m_el = std::get<1>(data); UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); Real M_n = rho * A * L / 2; m_el(0) = m_el(1) = M_n; } this->getDOFManager().assembleElementalArrayLocalArray( m_all_el, M, _segment_2, ghost_type); this->getDOFManager().assembleToLumpedMatrix("disp", M, "M"); } void assembleMass() { SparseMatrix & M = this->getDOFManager().getMatrix("M"); M.zero(); Array m_all_el(this->nb_elements, 4); Matrix m(2, 2); m(0, 0) = m(1, 1) = 2; m(0, 1) = m(1, 0) = 1; // under integrated // m(0, 0) = m(1, 1) = 3./2.; // m(0, 1) = m(1, 0) = 3./2.; // lumping the mass matrix // m(0, 0) += m(0, 1); // m(1, 1) += m(1, 0); // m(0, 1) = m(1, 0) = 0; for (auto && data : zip(make_view(this->mesh.getConnectivity(_segment_2), 2), make_view(m_all_el, 2, 2))) { const auto & conn = std::get<0>(data); auto & m_el = std::get<1>(data); UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); m_el = m; m_el *= rho * A * L / 6.; } this->getDOFManager().assembleElementalMatricesToMatrix( "M", "disp", m_all_el, _segment_2); is_mass_assembled = true; } - MatrixType getMatrixType(const ID &) override { return _symmetric; } + MatrixType getMatrixType(const ID &) const override { return _symmetric; } void assembleMatrix(const ID & matrix_id) override { if (matrix_id == "K") { if (not is_stiffness_assembled) this->assembleStiffness(); } else if (matrix_id == "M") { if (not is_mass_assembled) this->assembleMass(); } else if (matrix_id == "C") { // pass, no damping matrix } else { AKANTU_EXCEPTION("This solver does not know what to do with a matrix " << matrix_id); } } void assembleLumpedMatrix(const ID & matrix_id) override { if (matrix_id == "M") { if (not is_lumped_mass_assembled) this->assembleLumpedMass(); } else { AKANTU_EXCEPTION("This solver does not know what to do with a matrix " << matrix_id); } } void assembleStiffness() { SparseMatrix & K = this->getDOFManager().getMatrix("K"); K.zero(); Matrix k(2, 2); k(0, 0) = k(1, 1) = 1; k(0, 1) = k(1, 0) = -1; Array k_all_el(this->nb_elements, 4); auto k_it = k_all_el.begin(2, 2); auto cit = this->mesh.getConnectivity(_segment_2).begin(2); auto cend = this->mesh.getConnectivity(_segment_2).end(2); for (; cit != cend; ++cit, ++k_it) { const auto & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); auto & k_el = *k_it; k_el = k; k_el *= E * A / L; } this->getDOFManager().assembleElementalMatricesToMatrix( "K", "disp", k_all_el, _segment_2); is_stiffness_assembled = true; } void assembleResidual() override { this->getDOFManager().assembleToResidual("disp", forces); internal_forces.zero(); this->assembleResidualInternal(_not_ghost); this->synchronize(SynchronizationTag::_user_1); this->getDOFManager().assembleToResidual("disp", internal_forces, -1.); } void assembleResidualInternal(GhostType ghost_type) { Array forces_internal_el( this->mesh.getNbElement(_segment_2, ghost_type), 2); auto cit = this->mesh.getConnectivity(_segment_2, ghost_type).begin(2); auto cend = this->mesh.getConnectivity(_segment_2, ghost_type).end(2); auto f_it = forces_internal_el.begin(2); auto strain_it = this->strains.begin(); auto stress_it = this->stresses.begin(); auto L_it = this->initial_lengths.begin(); for (; cit != cend; ++cit, ++f_it, ++strain_it, ++stress_it, ++L_it) { const auto & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real u1 = this->displacement(n1, _x); Real u2 = this->displacement(n2, _x); *strain_it = (u2 - u1) / *L_it; *stress_it = E * *strain_it; Real f_n = A * *stress_it; Vector & f = *f_it; f(0) = -f_n; f(1) = f_n; } this->getDOFManager().assembleElementalArrayLocalArray( forces_internal_el, internal_forces, _segment_2, ghost_type); } Real getPotentialEnergy() { Real res = 0; if (not lumped) { res = this->mulVectMatVect(this->displacement, "K", this->displacement); } else { auto strain_it = this->strains.begin(); auto stress_it = this->stresses.begin(); auto strain_end = this->strains.end(); auto L_it = this->initial_lengths.begin(); for (; strain_it != strain_end; ++strain_it, ++stress_it, ++L_it) { res += *strain_it * *stress_it * A * *L_it; } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); } return res / 2.; } Real getKineticEnergy() { Real res = 0; if (not lumped) { res = this->mulVectMatVect(this->velocity, "M", this->velocity); } else { Array & m = dynamic_cast( this->getDOFManager().getLumpedMatrix("M")); auto it = velocity.begin(); auto end = velocity.end(); auto m_it = m.begin(); for (UInt node = 0; it != end; ++it, ++m_it, ++node) { if (mesh.isLocalOrMasterNode(node)) res += *m_it * *it * *it; } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); } return res / 2.; } Real getExternalWorkIncrement() { Real res = 0; auto it = velocity.begin(); auto end = velocity.end(); auto if_it = internal_forces.begin(); auto ef_it = forces.begin(); auto b_it = blocked.begin(); for (UInt node = 0; it != end; ++it, ++if_it, ++ef_it, ++b_it, ++node) { if (mesh.isLocalOrMasterNode(node)) res += (*b_it ? -*if_it : *ef_it) * *it; } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); return res * this->getTimeStep(); } Real mulVectMatVect(const Array & x, const ID & A_id, const Array & y) { Array Ay(nb_dofs); this->getDOFManager().assembleMatMulVectToArray("disp", A_id, y, Ay); Real res = 0.; for (auto && data : zip(arange(nb_dofs), make_view(Ay), make_view(x))) { res += std::get<2>(data) * std::get<1>(data) * mesh.isLocalOrMasterNode(std::get<0>(data)); } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); return res; } /* ------------------------------------------------------------------------ */ UInt getNbData(const Array & elements, const SynchronizationTag &) const override { return elements.size() * sizeof(Real); } void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_user_1) { for (const auto & el : elements) { buffer << this->stresses(el.element); } } } void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override { if (tag == SynchronizationTag::_user_1) { auto cit = this->mesh.getConnectivity(_segment_2, _ghost).begin(2); for (const auto & el : elements) { Real stress; buffer >> stress; Real f = A * stress; Vector conn = cit[el.element]; this->internal_forces(conn(0), _x) += -f; this->internal_forces(conn(1), _x) += f; } } } const Mesh & getMesh() const { return mesh; } UInt getSpatialDimension() const { return 1; } auto & getBlockedDOFs() { return blocked; } private: UInt nb_dofs; UInt nb_elements; bool lumped; bool is_stiffness_assembled{false}; bool is_mass_assembled{false}; bool is_lumped_mass_assembled{false}; public: Real E, A, rho; Mesh & mesh; Array displacement; Array velocity; Array acceleration; Array blocked; Array forces; Array internal_forces; Array stresses; Array strains; Array initial_lengths; }; #endif /* AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH_ */ } // namespace akantu diff --git a/test/test_model/test_common/test_non_local_toolbox/my_model.hh b/test/test_model/test_common/test_non_local_toolbox/my_model.hh index 4dba337a2..9c0fa5e3e 100644 --- a/test/test_model/test_common/test_non_local_toolbox/my_model.hh +++ b/test/test_model/test_common/test_non_local_toolbox/my_model.hh @@ -1,124 +1,126 @@ /** * @file my_model.hh * * @author Nicolas Richart * * @date creation: Mon Sep 11 2017 * @date last modification: Fri Jun 26 2020 * * @brief A dummy model for tests purposes * * * @section LICENSE * * Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "integrator_gauss.hh" #include "model.hh" #include "non_local_manager.hh" #include "non_local_manager_callback.hh" #include "non_local_neighborhood_base.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; class MyModel : public Model, public NonLocalManagerCallback { using MyFEEngineType = FEEngineTemplate; public: MyModel(Mesh & mesh, UInt spatial_dimension) : Model(mesh, ModelType::_model, spatial_dimension), manager(*this, *this) { registerFEEngineObject("FEEngine", mesh, spatial_dimension); manager.registerNeighborhood("test_region", "test_region"); getFEEngine().initShapeFunctions(); manager.initialize(); } void initModel() override {} - MatrixType getMatrixType(const ID &) override { return _mt_not_defined; } + MatrixType getMatrixType(const ID &) const override { + return _mt_not_defined; + } std::tuple getDefaultSolverID(const AnalysisMethod & /*method*/) override { return std::make_tuple("test", TimeStepSolverType::_static); } void assembleMatrix(const ID &) override {} void assembleLumpedMatrix(const ID &) override {} void assembleResidual() override {} void onNodesAdded(const Array &, const NewNodesEvent &) override {} void onNodesRemoved(const Array &, const Array &, const RemovedNodesEvent &) override {} void onElementsAdded(const Array &, const NewElementsEvent &) override {} void onElementsRemoved(const Array &, const ElementTypeMapArray &, const RemovedElementsEvent &) override {} void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override {} void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override { ElementTypeMapArray quadrature_points_coordinates( "quadrature_points_coordinates_tmp_nl", this->id); quadrature_points_coordinates.initialize(this->getFEEngine(), _nb_component = spatial_dimension, _ghost_type = ghost_type); IntegrationPoint q; q.ghost_type = ghost_type; q.global_num = 0; auto & neighborhood = manager.getNeighborhood("test_region"); for (auto & type : quadrature_points_coordinates.elementTypes( spatial_dimension, ghost_type)) { q.type = type; auto & quads = quadrature_points_coordinates(type, ghost_type); this->getFEEngine().computeIntegrationPointsCoordinates(quads, type, ghost_type); auto quad_it = quads.begin(quads.getNbComponent()); auto quad_end = quads.end(quads.getNbComponent()); q.num_point = 0; for (; quad_it != quad_end; ++quad_it) { neighborhood.insertIntegrationPoint(q, *quad_it); ++q.num_point; ++q.global_num; } } } void computeNonLocalStresses(GhostType) override {} void updateLocalInternal(ElementTypeMapReal &, GhostType, ElementKind) override {} void updateNonLocalInternal(ElementTypeMapReal &, GhostType, ElementKind) override {} const auto & getNonLocalManager() const { return manager; } private: NonLocalManager manager; }; diff --git a/test/test_solver/test_sparse_solver_mumps.cc b/test/test_solver/test_sparse_solver_mumps.cc index a96055b2b..ed881567c 100644 --- a/test/test_solver/test_sparse_solver_mumps.cc +++ b/test/test_solver/test_sparse_solver_mumps.cc @@ -1,169 +1,169 @@ /** * @file test_sparse_solver_mumps.cc * * @author Nicolas Richart * * @date creation: Fri May 19 2017 * @date last modification: Sun Dec 30 2018 * * @brief test the matrix vector product in parallel * * * @section LICENSE * * Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "dof_synchronizer.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "mesh_partition_scotch.hh" #include "sparse_matrix_aij.hh" #include "sparse_solver_mumps.hh" #include "terms_to_assemble.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); const UInt spatial_dimension = 1; const UInt nb_global_dof = 11; const auto & comm = Communicator::getStaticCommunicator(); Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); Mesh mesh(spatial_dimension); if (prank == 0) { genMesh(mesh, nb_global_dof); RandomGenerator::seed(1496137735); } else { RandomGenerator::seed(2992275470); } mesh.distribute(); UInt node = 0; for (auto pos : mesh.getNodes()) { std::cout << prank << " " << node << " pos: " << pos << " [" << mesh.getNodeGlobalId(node) << "] " << mesh.getNodeFlag(node) << std::endl; ++node; } UInt nb_nodes = mesh.getNbNodes(); DOFManagerDefault dof_manager(mesh, "test_dof_manager"); Array x(nb_nodes); dof_manager.registerDOFs("x", x, _dst_nodal); const auto & local_equation_number = dof_manager.getLocalEquationsNumbers("x"); auto & A = dof_manager.getNewMatrix("A", _symmetric); Array b(nb_nodes); - TermsToAssemble terms; + TermsToAssemble terms("x", "x"); for (UInt i = 0; i < nb_nodes; ++i) { if (dof_manager.isLocalOrMasterDOF(i)) { auto li = local_equation_number(i); auto gi = dof_manager.localToGlobalEquationNumber(li); terms(i, i) = 1. / (1. + gi); } } - dof_manager.assemblePreassembledMatrix("x", "x", "A", terms); + dof_manager.assemblePreassembledMatrix("A", terms); std::stringstream str; str << "Matrix_" << prank << ".mtx"; A.saveMatrix(str.str()); for (UInt n = 0; n < nb_nodes; ++n) { b(n) = 1.; } SparseSolverMumps solver(dof_manager, "A"); solver.solve(x, b); auto && check = [&](auto && xs) { debug::setDebugLevel(dblTest); std::cout << xs << std::endl; debug::setDebugLevel(dblWarning); UInt d = 1.; for (auto x : xs) { if (std::abs(x - d) / d > 1e-15) AKANTU_EXCEPTION("Error in the solution: " << x << " != " << d << " [" << (std::abs(x - d) / d) << "]."); ++d; } }; if (psize > 1) { auto & sync = dynamic_cast(dof_manager).getSynchronizer(); if (prank == 0) { Array x_gathered(dof_manager.getSystemSize()); sync.gather(x, x_gathered); check(x_gathered); } else { sync.gather(x); } } else { check(x); } finalize(); return 0; } /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes) { MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); Array & conn = mesh_accessor.getConnectivity(_segment_2); nodes.resize(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); } conn.resize(nb_nodes - 1); for (UInt n = 0; n < nb_nodes - 1; ++n) { conn(n, 0) = n; conn(n, 1) = n + 1; } mesh_accessor.makeReady(); }