diff --git a/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index 5c8cdcbe1..5b79a3878 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,813 +1,813 @@ /** * @file heat_transfer_model.cc * * @author Rui Wang * @author Srinivasa Babu Ramisetti * @author Guillaume Anciaux * @author Nicolas Richart * @author David Simon Kammer * * @date Sun May 01 19:14:43 2011 * * @brief Implementation of HeatTransferModel class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "heat_transfer_model.hh" #include "aka_math.hh" #include "aka_common.hh" #include "fem_template.hh" #include "mesh.hh" #include "static_communicator.hh" #include "parser.hh" #include "generalized_trapezoidal.hh" #ifdef AKANTU_USE_IOHELPER -# include "dumper_iohelper.hh" +# include "dumper_paraview.hh" # include "dumper_iohelper_tmpl_homogenizing_field.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ HeatTransferModel::HeatTransferModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), Dumpable(), integrator(new ForwardEuler()), temperature_gradient ("temperature_gradient", id), temperature_on_qpoints ("temperature_on_qpoints", id), conductivity_on_qpoints ("conductivity_on_qpoints", id), k_gradt_on_qpoints ("k_gradt_on_qpoints", id), int_bt_k_gT ("int_bt_k_gT", id), bt_k_gT ("bt_k_gT", id), conductivity(spatial_dimension, spatial_dimension), thermal_energy ("thermal_energy", id) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); std::stringstream sstr; sstr << id << ":fem"; registerFEMObject(sstr.str(), mesh,spatial_dimension); this->temperature= NULL; this->residual = NULL; this->boundary = NULL; this->conductivity_variation = 0.0; this->T_ref = 0; this->registerDumper("paraview_all", id, true); this->addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initModel() { getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; Synchronizer & synch_parallel = createParallelSynch(partition,data_accessor); synch_registry->registerSynchronizer(synch_parallel, _gst_htm_capacity); synch_registry->registerSynchronizer(synch_parallel, _gst_htm_temperature); synch_registry->registerSynchronizer(synch_parallel, _gst_htm_gradient_temperature); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initPBC() { AKANTU_DEBUG_IN(); Model::initPBC(); PBCSynchronizer * synch = new PBCSynchronizer(pbc_pair); synch_registry->registerSynchronizer(*synch, _gst_htm_capacity); synch_registry->registerSynchronizer(*synch, _gst_htm_temperature); changeLocalEquationNumberforPBC(pbc_pair,1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); std::stringstream sstr_temp; sstr_temp << id << ":temperature"; std::stringstream sstr_temp_rate; sstr_temp_rate << id << ":temperature_rate"; std::stringstream sstr_inc; sstr_inc << id << ":increment"; std::stringstream sstr_ext_flx; sstr_ext_flx << id << ":external_flux"; std::stringstream sstr_residual; sstr_residual << id << ":residual"; std::stringstream sstr_lump; sstr_lump << id << ":lumped"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; temperature = &(alloc(sstr_temp.str(), nb_nodes, 1, REAL_INIT_VALUE)); temperature_rate = &(alloc(sstr_temp_rate.str(), nb_nodes, 1, REAL_INIT_VALUE)); increment = &(alloc(sstr_inc.str(), nb_nodes, 1, REAL_INIT_VALUE)); external_heat_rate = &(alloc(sstr_ext_flx.str(), nb_nodes, 1, REAL_INIT_VALUE)); residual = &(alloc(sstr_residual.str(), nb_nodes, 1, REAL_INIT_VALUE)); capacity_lumped = &(alloc(sstr_lump.str(), nb_nodes, 1, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, 1, false)); Mesh::ConnectivityTypeList::const_iterator it; /* -------------------------------------------------------------------------- */ // byelementtype vectors getFEM().getMesh().initByElementTypeArray(temperature_on_qpoints, 1, spatial_dimension); getFEM().getMesh().initByElementTypeArray(temperature_gradient, spatial_dimension, spatial_dimension); getFEM().getMesh().initByElementTypeArray(conductivity_on_qpoints, spatial_dimension*spatial_dimension, spatial_dimension); getFEM().getMesh().initByElementTypeArray(k_gradt_on_qpoints, spatial_dimension, spatial_dimension); getFEM().getMesh().initByElementTypeArray(bt_k_gT, 1, spatial_dimension,true); getFEM().getMesh().initByElementTypeArray(int_bt_k_gT, 1, spatial_dimension,true); getFEM().getMesh().initByElementTypeArray(thermal_energy, 1, spatial_dimension); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(gt); for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_element = getFEM().getMesh().getNbElement(*it, gt); UInt nb_quad_points = this->getFEM().getNbQuadraturePoints(*it, gt) * nb_element; temperature_on_qpoints(*it, gt).resize(nb_quad_points); temperature_on_qpoints(*it, gt).clear(); temperature_gradient(*it, gt).resize(nb_quad_points); temperature_gradient(*it, gt).clear(); conductivity_on_qpoints(*it, gt).resize(nb_quad_points); conductivity_on_qpoints(*it, gt).clear(); k_gradt_on_qpoints(*it, gt).resize(nb_quad_points); k_gradt_on_qpoints(*it, gt).clear(); bt_k_gT(*it, gt).resize(nb_quad_points); bt_k_gT(*it, gt).clear(); int_bt_k_gT(*it, gt).resize(nb_element); int_bt_k_gT(*it, gt).clear(); thermal_energy(*it, gt).resize(nb_element); thermal_energy(*it, gt).clear(); } } /* -------------------------------------------------------------------------- */ dof_synchronizer = new DOFSynchronizer(getFEM().getMesh(),1); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ HeatTransferModel::~HeatTransferModel() { AKANTU_DEBUG_IN(); if(dof_synchronizer) delete dof_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); FEM & fem = getFEM(); const Mesh::ConnectivityTypeList & type_list = fem.getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_element = getFEM().getMesh().getNbElement(*it,ghost_type); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, ghost_type); Array rho_1 (nb_element * nb_quadrature_points,1, capacity * density); fem.assembleFieldLumped(rho_1,1,*capacity_lumped, dof_synchronizer->getLocalDOFEquationNumbers(), *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped() { AKANTU_DEBUG_IN(); capacity_lumped->clear(); assembleCapacityLumped(_not_ghost); assembleCapacityLumped(_ghost); getSynchronizerRegistry().synchronize(_gst_htm_capacity); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::updateResidual() { AKANTU_DEBUG_IN(); /// @f$ r = q_{ext} - q_{int} - C \dot T @f$ // start synchronization synch_registry->asynchronousSynchronize(_gst_htm_temperature); // finalize communications synch_registry->waitEndSynchronize(_gst_htm_temperature); //clear the array /// first @f$ r = q_{ext} @f$ // residual->clear(); residual->copy(*external_heat_rate); /// then @f$ r -= q_{int} @f$ // update the not ghost ones updateResidual(_not_ghost); // update for the received ghosts updateResidual(_ghost); /// finally @f$ r -= C \dot T @f$ // lumped C UInt nb_nodes = temperature_rate->getSize(); UInt nb_degree_of_freedom = temperature_rate->getNbComponent(); Real * capacity_val = capacity_lumped->values; Real * temp_rate_val = temperature_rate->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *res_val -= *capacity_val * *temp_rate_val; } boundary_val++; res_val++; capacity_val++; temp_rate_val++; } #ifndef AKANTU_NDEBUG getSynchronizerRegistry().synchronize(akantu::_gst_htm_gradient_temperature); #endif solveExplicitLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeConductivityOnQuadPoints(const GhostType & ghost_type) { const Mesh::ConnectivityTypeList & type_list = this->getFEM().getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; Array & temperature_interpolated = temperature_on_qpoints(*it, ghost_type); //compute the temperature on quadrature points this->getFEM().interpolateOnQuadraturePoints(*temperature, temperature_interpolated, 1 ,*it,ghost_type); Array::iterator< Matrix > C_it = conductivity_on_qpoints(*it, ghost_type).begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > C_end = conductivity_on_qpoints(*it, ghost_type).end(spatial_dimension, spatial_dimension); Array::iterator T_it = temperature_interpolated.begin(); for (;C_it != C_end; ++C_it, ++T_it) { Matrix & C = *C_it; Real & T = *T_it; C = conductivity; Matrix variation(spatial_dimension, spatial_dimension, conductivity_variation * (T - T_ref)); C += conductivity_variation; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeKgradT(const GhostType & ghost_type) { computeConductivityOnQuadPoints(ghost_type); const Mesh::ConnectivityTypeList & type_list = this->getFEM().getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { const ElementType & type = *it; if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; Array & gradient = temperature_gradient(*it, ghost_type); this->getFEM().gradientOnQuadraturePoints(*temperature, gradient, 1 ,*it,ghost_type); Array::iterator< Matrix > C_it = conductivity_on_qpoints(*it, ghost_type).begin(spatial_dimension, spatial_dimension); Array::iterator< Vector > BT_it = gradient.begin(spatial_dimension); Array::iterator< Vector > k_BT_it = k_gradt_on_qpoints(type, ghost_type).begin(spatial_dimension); Array::iterator< Vector > k_BT_end = k_gradt_on_qpoints(type, ghost_type).end(spatial_dimension); for (;k_BT_it != k_BT_end; ++k_BT_it, ++BT_it, ++C_it) { Vector & k_BT = *k_BT_it; Vector & BT = *BT_it; Matrix & C = *C_it; k_BT.mul(C, BT); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::updateResidual(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const Mesh::ConnectivityTypeList & type_list = this->getFEM().getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; Array & shapes_derivatives = const_cast &>(getFEM().getShapesDerivatives(*it,ghost_type)); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); // compute k \grad T computeKgradT(ghost_type); Array::iterator< Vector > k_BT_it = k_gradt_on_qpoints(*it,ghost_type).begin(spatial_dimension); Array::iterator< Matrix > B_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::iterator< Vector > Bt_k_BT_it = bt_k_gT(*it,ghost_type).begin(nb_nodes_per_element); Array::iterator< Vector > Bt_k_BT_end = bt_k_gT(*it,ghost_type).end(nb_nodes_per_element); for (;Bt_k_BT_it != Bt_k_BT_end; ++Bt_k_BT_it, ++B_it, ++k_BT_it) { Vector & k_BT = *k_BT_it; Vector & Bt_k_BT = *Bt_k_BT_it; Matrix & B = *B_it; Bt_k_BT.mul(B, k_BT); } this->getFEM().integrate(bt_k_gT(*it,ghost_type), int_bt_k_gT(*it,ghost_type), nb_nodes_per_element, *it,ghost_type); this->getFEM().assembleArray(int_bt_k_gT(*it,ghost_type), *residual, dof_synchronizer->getLocalDOFEquationNumbers(), 1, *it,ghost_type, empty_filter, -1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::solveExplicitLumped() { AKANTU_DEBUG_IN(); UInt nb_nodes = increment->getSize(); UInt nb_degree_of_freedom = increment->getNbComponent(); Real * capa_val = capacity_lumped->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; Real * inc = increment->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *inc = (*res_val / *capa_val); } res_val++; boundary_val++; inc++; capa_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::explicitPred() { AKANTU_DEBUG_IN(); integrator->integrationSchemePred(time_step, *temperature, *temperature_rate, *boundary); UInt nb_nodes = temperature->getSize(); UInt nb_degree_of_freedom = temperature->getNbComponent(); Real * temp = temperature->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n, ++temp) if(*temp < 0.) *temp = 0.; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrTempRate(time_step, *temperature, *temperature_rate, *boundary, *increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real el_size; Real min_el_size = std::numeric_limits::max(); Real conductivitymax = conductivity(0, 0); //get the biggest parameter from k11 until k33// for(UInt i = 0; i < spatial_dimension; i++) for(UInt j = 0; j < spatial_dimension; j++) conductivitymax = std::max(conductivity(i, j), conductivitymax); const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_nodes_per_element = getFEM().getMesh().getNbNodesPerElement(*it); Array coord(0, nb_nodes_per_element*spatial_dimension); FEM::extractNodalToElementField(getFEM().getMesh(), getFEM().getMesh().getNodes(), coord, *it, _not_ghost); Array::iterator< Matrix > el_coord = coord.begin(spatial_dimension, nb_nodes_per_element); UInt nb_element = getFEM().getMesh().getNbElement(*it); for (UInt el = 0; el < nb_element; ++el, ++el_coord) { el_size = getFEM().getElementInradius(*el_coord, *it); min_el_size = std::min(min_el_size, el_size); } AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size << " and the max conductivity is : " << conductivitymax); } Real min_dt = 2 * min_el_size * min_el_size * density * capacity / conductivitymax; StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::readMaterials(const std::string & filename) { Parser parser; parser.open(filename); std::string opt_param; std::string mat_type = parser.getNextSection("heat", opt_param); if (mat_type != "") { parser.readSection(*this); } else AKANTU_DEBUG_ERROR("did not find any section with material info " <<"for heat conduction"); } /* -------------------------------------------------------------------------- */ bool HeatTransferModel::setParam(const std::string & key, const std::string & value) { std::stringstream str(value); if (key == "conductivity") { for(UInt i = 0; i < 3; i++) for(UInt j = 0; j < 3; j++) { if (i< spatial_dimension && j < spatial_dimension){ str >> conductivity(i, j); AKANTU_DEBUG_INFO("Conductivity(" << i << "," << j << ") = " << conductivity[i*spatial_dimension+j]); } else { Real tmp; str >> tmp; } } } else if (key == "conductivity_variation") { str >> conductivity_variation; } else if (key == "temperature_reference") { str >> T_ref; } else if (key == "capacity"){ str >> capacity; AKANTU_DEBUG_INFO("The capacity of the material is:" << capacity); } else if (key == "density"){ str >> density; AKANTU_DEBUG_INFO("The density of the material is:" << density); } else { return false; } return true; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initFull(const std::string & material_file){ readMaterials(material_file); //model initialization initModel(); //initialize the vectors initArrays(); temperature->clear(); temperature_rate->clear(); external_heat_rate->clear(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initFEMBoundary(bool create_surface) { if(create_surface) MeshUtils::buildFacets(getFEM().getMesh()); FEM & fem_boundary = getFEMBoundary(); fem_boundary.initShapeFunctions(); fem_boundary.computeNormalsOnControlPoints(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::computeThermalEnergyByNode() { AKANTU_DEBUG_IN(); Real ethermal = 0.; Array::iterator< Vector > heat_rate_it = residual->begin(residual->getNbComponent()); Array::iterator< Vector > heat_rate_end = residual->end(residual->getNbComponent()); UInt n = 0; for(;heat_rate_it != heat_rate_end; ++heat_rate_it, ++n) { Real heat = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; Vector & heat_rate = *heat_rate_it; for (UInt i = 0; i < heat_rate.size(); ++i) { if (count_node) heat += heat_rate[i] * time_step; } ethermal += heat; } StaticCommunicator::getStaticCommunicator().allReduce(ðermal, 1, _so_sum); AKANTU_DEBUG_OUT(); return ethermal; } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::getThermalEnergy(iterator Eth, Array::const_iterator T_it, Array::const_iterator T_end) const { for(;T_it != T_end; ++T_it, ++Eth) { *Eth = capacity * density * *T_it; } } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(type); Vector Eth_on_quarature_points(nb_quadrature_points); Array::iterator T_it = this->temperature_on_qpoints(type).begin(); T_it += index * nb_quadrature_points; Array::iterator T_end = T_it + nb_quadrature_points; getThermalEnergy(Eth_on_quarature_points.storage(), T_it, T_end); return getFEM().integrate(Eth_on_quarature_points, type, index); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy() { Real Eth = 0; Mesh & mesh = getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { UInt nb_element = getFEM().getMesh().getNbElement(*it, _not_ghost); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, _not_ghost); Array Eth_per_quad(nb_element * nb_quadrature_points, 1); Array::iterator T_it = this->temperature_on_qpoints(*it).begin(); Array::iterator T_end = this->temperature_on_qpoints(*it).end(); getThermalEnergy(Eth_per_quad.begin(), T_it, T_end); Eth += getFEM().integrate(Eth, *it); } return Eth; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id) { AKANTU_DEBUG_IN(); Real energy = 0; if("thermal") energy = getThermalEnergy(); // reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real energy = 0.; if("thermal") energy = getThermalEnergy(type, index); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #define IF_ADD_NODAL_FIELD(dumper_name, field, type, pad) \ if(field_id == BOOST_PP_STRINGIZE(field)) { \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field); \ if(pad) f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), f); \ } else #define IF_ADD_ELEM_FIELD(dumper_name, field, type, pad) \ if(field_id == BOOST_PP_STRINGIZE(field)) { \ typedef DumperIOHelper::HomogenizedField Field; \ DumperIOHelper::Field * f = \ new Field(getFEM(), \ field, \ spatial_dimension, \ spatial_dimension, \ _not_ghost, \ _ek_regular); \ if(pad) f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), f); \ } else #endif void HeatTransferModel::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER IF_ADD_NODAL_FIELD(dumper_name, temperature , Real, false) IF_ADD_NODAL_FIELD(dumper_name, temperature_rate , Real, false) IF_ADD_NODAL_FIELD(dumper_name, external_heat_rate, Real, false) IF_ADD_NODAL_FIELD(dumper_name, residual , Real, false) IF_ADD_NODAL_FIELD(dumper_name, capacity_lumped , Real, false) IF_ADD_NODAL_FIELD(dumper_name, boundary , bool, false) IF_ADD_ELEM_FIELD (dumper_name, temperature_gradient, Real, false) if(field_id == "partitions" ) { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField<>(mesh, spatial_dimension, _not_ghost, _ek_regular)); } else { AKANTU_DEBUG_ERROR("Field " << field_id << " does not exists in the model " << id); } #endif } /* -------------------------------------------------------------------------- */ void HeatTransferModel::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER IF_ADD_ELEM_FIELD (dumper_name, temperature_gradient, Real, false) { AKANTU_DEBUG_ERROR("Field " << field_id << " does not exists in the model " << id << " or is not dumpable as a vector"); } #endif } /* -------------------------------------------------------------------------- */ void HeatTransferModel::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER AKANTU_DEBUG_ERROR("Field " << field_id << " does not exists in the model " << id << " or is not dumpable as a Tensor"); #endif } __END_AKANTU__ diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index dd4185463..a691bd816 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,657 +1,661 @@ /** * @file structural_mechanics_model.cc * @author Fabian Barras * @date Thu May 5 15:52:38 2011 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 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 "aka_math.hh" #include "integration_scheme_2nd_order.hh" #include "static_communicator.hh" #include "sparse_matrix.hh" #include "solver.hh" #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif + +#ifdef AKANTU_USE_IOHELPER +# include "dumper_paraview.hh" +#endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), Dumpable(), stress("stress", id, memory_id), element_material("element_material", id, memory_id), stiffness_matrix(NULL), solver(NULL), rotation_matrix("rotation_matices", id, memory_id) { AKANTU_DEBUG_IN(); registerFEMObject("StructuralMechanicsFEM", mesh, spatial_dimension); this->displacement_rotation = NULL; this->force_momentum = NULL; this->residual = NULL; this->boundary = NULL; this->increment = NULL; if(spatial_dimension == 2) nb_degree_of_freedom = 3; else if (spatial_dimension == 3) nb_degree_of_freedom = 6; else { AKANTU_DEBUG_TO_IMPLEMENT(); } this->registerDumper("paraview_all", id, true); this->addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() { AKANTU_DEBUG_IN(); if(solver) delete solver; if(stiffness_matrix) delete stiffness_matrix; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFull(std::string material) { initModel(); initArrays(); initSolver(); displacement_rotation->clear(); force_momentum ->clear(); residual ->clear(); boundary ->clear(); increment ->clear(); Mesh::type_iterator it = getFEM().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEM().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { computeRotationMatrix(*it); } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; std::stringstream sstr_incr; sstr_incr << id << ":increment"; displacement_rotation = &(alloc(sstr_disp.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); force_momentum = &(alloc(sstr_forc.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, nb_degree_of_freedom, false)); increment = &(alloc(sstr_incr.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); Mesh::type_iterator it = getFEM().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEM().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { UInt nb_element = getFEM().getMesh().getNbElement(*it); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it); element_material.alloc(nb_element, 1, *it, _not_ghost); set_ID.alloc(nb_element, 1, *it, _not_ghost); UInt size = getTangentStiffnessVoigtSize(*it); stress.alloc(nb_element * nb_quadrature_points, size , *it, _not_ghost); } dof_synchronizer = new DOFSynchronizer(getFEM().getMesh(), nb_degree_of_freedom); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { getFEM().initShapeFunctions(_not_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { AKANTU_DEBUG_IN(); const Mesh & mesh = getFEM().getMesh(); #if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_node = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(nb_global_node * nb_degree_of_freedom, _symmetric, nb_degree_of_freedom, sstr.str(), memory_id); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS solver->initialize(options); #endif //AKANTU_HAS_SOLVER AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt StructuralMechanicsModel::getTangentStiffnessVoigtSize(const ElementType & type) { UInt size; #define GET_TANGENT_STIFFNESS_VOIGT_SIZE(type) \ size = getTangentStiffnessVoigtSize(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_TANGENT_STIFFNESS_VOIGT_SIZE); #undef GET_TANGENT_STIFFNESS_VOIGT_SIZE return size; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); stiffness_matrix->clear(); Mesh::type_iterator it = getFEM().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEM().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { ElementType type = *it; #define ASSEMBLE_STIFFNESS_MATRIX(type) \ assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_2>(Array & rotations){ ElementType type = _bernoulli_beam_2; Mesh & mesh = getFEM().getMesh(); UInt nb_element = mesh.getNbElement(type); Array::iterator< Vector > connec_it = mesh.getConnectivity(type).begin(2); Array::iterator< Vector > nodes_it = mesh.getNodes().begin(spatial_dimension); Array::iterator< Matrix > R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); for (UInt e = 0; e < nb_element; ++e, ++R_it, ++connec_it) { Matrix & R = *R_it; Vector & connec = *connec_it; Vector x2 = nodes_it[connec(1)]; // X2 Vector x1 = nodes_it[connec(0)]; // X1 Real le = x1.distance(x2); Real c = (x2(0) - x1(0)) / le; Real s = (x2(1) - x1(1)) / le; /// Definition of the rotation matrix R(0,0) = c; R(0,1) = s; R(0,2) = 0.; R(1,0) = -s; R(1,1) = c; R(1,2) = 0.; R(2,0) = 0.; R(2,1) = 0.; R(2,2) = 1.; } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_3>(Array & rotations){ ElementType type = _bernoulli_beam_3; Mesh & mesh = getFEM().getMesh(); UInt nb_element = mesh.getNbElement(type); Array::iterator< Vector > n_it = mesh.getNormals(type).begin(spatial_dimension); Array::iterator< Vector > connec_it = mesh.getConnectivity(type).begin(2); Array::iterator< Vector > nodes_it = mesh.getNodes().begin(spatial_dimension); Matrix Pe (spatial_dimension, spatial_dimension); Matrix Pg (spatial_dimension, spatial_dimension); Matrix inv_Pg(spatial_dimension, spatial_dimension); Vector x_n(spatial_dimension); // x vect n Array::iterator< Matrix > R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); for (UInt e=0 ; e < nb_element; ++e, ++n_it, ++connec_it, ++R_it) { Vector & n = *n_it; Matrix & R = *R_it; Vector & connec = *connec_it; Vector x = nodes_it[connec(1)]; // X2 Vector y = nodes_it[connec(0)]; // X1 Real l = x.distance(y); x -= y; // X2 - X1 x_n.crossProduct(x, n); Pe.eye(); Pe(0, 0) *= l; Pe(1, 1) *= -l; Pg(0,0) = x(0); Pg(0,1) = x_n(0); Pg(0,2) = n(0); Pg(1,0) = x(1); Pg(1,1) = x_n(1); Pg(1,2) = n(1); Pg(2,0) = x(2); Pg(2,1) = x_n(2); Pg(2,2) = n(2); inv_Pg.inverse(Pg); Pe *= inv_Pg; for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { R(i, j) = Pe(i, j); R(i + spatial_dimension,j + spatial_dimension) = Pe(i, j); } } } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeRotationMatrix(const ElementType & type) { Mesh & mesh = getFEM().getMesh(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); if(!rotation_matrix.exists(type)) { rotation_matrix.alloc(nb_element, nb_degree_of_freedom*nb_nodes_per_element * nb_degree_of_freedom*nb_nodes_per_element, type); } else { rotation_matrix(type).resize(nb_element); } rotation_matrix(type).clear(); Arrayrotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom); rotations.clear(); #define COMPUTE_ROTATION_MATRIX(type) \ computeRotationMatrix(rotations); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX); #undef COMPUTE_ROTATION_MATRIX Array::iterator< Matrix > R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); Array::iterator< Matrix > T_it = rotation_matrix(type).begin(nb_degree_of_freedom*nb_nodes_per_element, nb_degree_of_freedom*nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) { Matrix & T = *T_it; Matrix & R = *R_it; T.clear(); for (UInt k = 0; k < nb_nodes_per_element; ++k){ for (UInt i = 0; i < nb_degree_of_freedom; ++i) for (UInt j = 0; j < nb_degree_of_freedom; ++j) T(k*nb_degree_of_freedom + i, k*nb_degree_of_freedom + j) = R(i, j); } } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); Mesh::type_iterator it = getFEM().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEM().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { ElementType type = *it; #define COMPUTE_STRESS_ON_QUAD(type) \ computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::updateResidual() { AKANTU_DEBUG_IN(); residual->copy(*force_momentum); Array ku(*displacement_rotation, true); ku *= *stiffness_matrix; *residual -= ku; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::solve() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving an implicit step."); UInt nb_nodes = displacement_rotation->getSize(); /// todo residual = force - Kxr * d_bloq jacobian_matrix->copyContent(*stiffness_matrix); jacobian_matrix->applyBoundary(*boundary); solver->setRHS(*residual); solver->solve(*increment); Real * increment_val = increment->values; Real * displacement_val = displacement_rotation->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *displacement_val += *increment_val; } displacement_val++; boundary_val++; increment_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance) { Real error; bool tmp = testConvergenceIncrement(tolerance, error); AKANTU_DEBUG_INFO("Norm of increment : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error) { AKANTU_DEBUG_IN(); Mesh & mesh= getFEM().getMesh(); UInt nb_nodes = displacement_rotation->getSize(); UInt nb_degree_of_freedom = displacement_rotation->getNbComponent(); Real norm = 0; Real * increment_val = increment->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*boundary_val) && is_local_node) { norm += *increment_val * *increment_val; } boundary_val++; increment_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); error = sqrt(norm); AKANTU_DEBUG_ASSERT(!isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (error < tolerance); } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_2>(Array & tangent_moduli) { UInt nb_element = getFEM().getMesh().getNbElement(_bernoulli_beam_2); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(_bernoulli_beam_2); UInt tangent_size = 2; Array::iterator< Matrix > D_it = tangent_moduli.begin(tangent_size, tangent_size); Array & el_mat = element_material(_bernoulli_beam_2, _not_ghost); for (UInt e = 0; e < nb_element; ++e) { UInt mat = el_mat(e); Real E = materials[mat].E; Real A = materials[mat].A; Real I = materials[mat].I; for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) { Matrix & D = *D_it; D(0,0) = E * A; D(1,1) = E * I; } } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_3>(Array & tangent_moduli) { UInt nb_element = getFEM().getMesh().getNbElement(_bernoulli_beam_3); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(_bernoulli_beam_3); UInt tangent_size = 4; Array::iterator< Matrix > D_it = tangent_moduli.begin(tangent_size, tangent_size); for (UInt e = 0; e < nb_element; ++e) { UInt mat = element_material(_bernoulli_beam_3, _not_ghost)(e); Real E = materials[mat].E; Real A = materials[mat].A; Real Iz = materials[mat].Iz; Real Iy = materials[mat].Iy; Real GJ = materials[mat].GJ; for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) { Matrix & D = *D_it; D(0,0) = E * A; D(1,1) = E * Iz; D(2,2) = E * Iy; D(3,3) = GJ; } } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_2>(Array & b, bool local) { UInt nb_element = getFEM().getMesh().getNbElement(_bernoulli_beam_2); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_2); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(_bernoulli_beam_2); MyFEMType & fem = getFEMClass(); Array::const_iterator< Vector > shape_Np = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 0).begin(nb_nodes_per_element); Array::const_iterator< Vector > shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 1).begin(nb_nodes_per_element); Array::const_iterator< Vector > shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 2).begin(nb_nodes_per_element); UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_2>(); UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom; b.clear(); Array::iterator< Matrix > B_it = b.begin(tangent_size, bt_d_b_size); for (UInt e = 0; e < nb_element; ++e) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Matrix & B = *B_it; const Vector & Np = *shape_Np; const Vector & Lpp = *shape_Lpp; const Vector & Mpp = *shape_Mpp; B(0,0) = Np(0); B(0,3) = Np(1); B(1,1) = Mpp(0); B(1,2) = Lpp(0); B(1,4) = Mpp(1); B(1,5) = Lpp(1); ++B_it; ++shape_Np; ++shape_Mpp; ++shape_Lpp; } // ++R_it; } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_3>(Array & b, bool local) { MyFEMType & fem = getFEMClass(); UInt nb_element = getFEM().getMesh().getNbElement(_bernoulli_beam_3); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_3); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(_bernoulli_beam_3); Array::const_iterator< Vector > shape_Np = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 0).begin(nb_nodes_per_element); Array::const_iterator< Vector > shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 1).begin(nb_nodes_per_element); Array::const_iterator< Vector > shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 2).begin(nb_nodes_per_element); UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_3>(); UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom; b.clear(); Array::iterator< Matrix > B_it = b.begin(tangent_size, bt_d_b_size); for (UInt e = 0; e < nb_element; ++e) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Matrix & B = *B_it; const Vector & Np = *shape_Np; const Vector & Lpp = *shape_Lpp; const Vector & Mpp = *shape_Mpp; B(0,0) = Np(0); B(0,6) = Np(1); B(1,1) = Mpp(0); B(1,5) = Lpp(0); B(1,7) = Mpp(1); B(1,11) = Lpp(1); B(2,2) = Mpp(0); B(2,4) = -Lpp(0); B(2,8) = Mpp(1); B(2,10) = -Lpp(1); B(3,3) = Np(0); B(3,9) = Np(1); ++B_it; ++shape_Np; ++shape_Mpp; ++shape_Lpp; } } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(dumper_name, id, field, type, n, stride) \ internalAddDumpFieldToDumper(dumper_name, id, \ new DumperIOHelper::NodalField(*field, n, stride)) UInt n; if(spatial_dimension == 2) { n = 2; } else n = 3; if(field_id == "displacement" ) { ADD_FIELD(dumper_name, "displacement", displacement_rotation, Real, n, 0); } else if(field_id == "rotation") { ADD_FIELD(dumper_name, "rotation", displacement_rotation, Real, nb_degree_of_freedom - n, n); } else if(field_id == "force" ) { ADD_FIELD(dumper_name, "force", force_momentum, Real, n, 0); } else if(field_id == "momentum") { ADD_FIELD(dumper_name, "momentum", force_momentum, Real, nb_degree_of_freedom - n, n); } else if(field_id == "residual") { ADD_FIELD(dumper_name, "residual", residual, Real, nb_degree_of_freedom, 0); } else if(field_id == "boundary") { ADD_FIELD(dumper_name, "boundary", boundary, bool, nb_degree_of_freedom, 0); } else if(field_id == "element_index_by_material") { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_material, spatial_dimension, _not_ghost, _ek_regular)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(dumper_name, id, field, type, n, stride) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field, n, stride); \ f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, id, f) UInt n; if(spatial_dimension == 2) { n = 2; } else n = 3; if(field_id == "displacement" ) { ADD_FIELD(dumper_name, "displacement", displacement_rotation, Real, n, 0); } else if(field_id == "force" ) { ADD_FIELD(dumper_name, "force", force_momentum, Real, n, 0); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::addDumpFieldTensorToDumper(__attribute__((unused)) const std::string & dumper_name, __attribute__((unused)) const std::string & field_id) { } __END_AKANTU__