diff --git a/src/fe_engine/shape_lagrange_inline_impl.cc b/src/fe_engine/shape_lagrange_inline_impl.cc index 6e9bbb1f6..a4b2cc5ea 100644 --- a/src/fe_engine/shape_lagrange_inline_impl.cc +++ b/src/fe_engine/shape_lagrange_inline_impl.cc @@ -1,527 +1,528 @@ /** * @file shape_lagrange_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Wed Oct 27 2010 * @date last modification: Thu Oct 15 2015 * * @brief ShapeLagrange inline implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" #include "aka_voigthelper.hh" #include "fe_engine.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__ #define __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ #define INIT_SHAPE_FUNCTIONS(type) \ setIntegrationPointsByType(integration_points, ghost_type); \ precomputeShapesOnIntegrationPoints(nodes, ghost_type); \ if (ElementClass::getNaturalSpaceDimension() == \ mesh.getSpatialDimension() || \ kind != _ek_regular) \ precomputeShapeDerivativesOnIntegrationPoints(nodes, ghost_type); template inline void ShapeLagrange::initShapeFunctions( const Array & nodes, const Matrix & integration_points, const ElementType & type, const GhostType & ghost_type) { AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS); } #undef INIT_SHAPE_FUNCTIONS /* -------------------------------------------------------------------------- */ template template inline void ShapeLagrange::computeShapeDerivativesOnCPointsByElement( const Matrix & node_coords, const Matrix & natural_coords, Tensor3 & shapesd) const { AKANTU_DEBUG_IN(); // compute dnds Tensor3 dnds(node_coords.rows(), node_coords.cols(), natural_coords.cols()); ElementClass::computeDNDS(natural_coords, dnds); // compute jacobian Tensor3 J(node_coords.rows(), natural_coords.rows(), natural_coords.cols()); ElementClass::computeJMat(dnds, node_coords, J); // compute dndx ElementClass::computeShapeDerivatives(J, dnds, shapesd); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::inverseMap(const Vector & real_coords, UInt elem, Vector & natural_coords, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(), elem_val + elem * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); ElementClass::inverseMap(real_coords, nodes_coord, natural_coords); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template bool ShapeLagrange::contains(const Vector & real_coords, UInt elem, const GhostType & ghost_type) const { UInt spatial_dimension = mesh.getSpatialDimension(); Vector natural_coords(spatial_dimension); inverseMap(real_coords, elem, natural_coords, ghost_type); return ElementClass::contains(natural_coords); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::interpolate(const Vector & real_coords, UInt elem, const Matrix & nodal_values, Vector & interpolated, const GhostType & ghost_type) const { UInt nb_shapes = ElementClass::getShapeSize(); Vector shapes(nb_shapes); computeShapes(real_coords, elem, shapes, ghost_type); ElementClass::interpolate(nodal_values, shapes, interpolated); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); Vector natural_coords(spatial_dimension); inverseMap(real_coords, elem, natural_coords, ghost_type); ElementClass::computeShapes(natural_coords, shapes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeShapeDerivatives( const Matrix & real_coords, UInt elem, Tensor3 & shapesd, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_points = real_coords.cols(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == shapesd.size(0) && nb_nodes_per_element == shapesd.size(1), "Shape size doesn't match"); AKANTU_DEBUG_ASSERT(nb_points == shapesd.size(2), "Number of points doesn't match shapes size"); Matrix natural_coords(spatial_dimension, nb_points); // Creates the matrix of natural coordinates for (UInt i = 0; i < nb_points; i++) { Vector real_point = real_coords(i); Vector natural_point = natural_coords(i); inverseMap(real_point, elem, natural_point, ghost_type); } UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(), elem_val + elem * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); computeShapeDerivativesOnCPointsByElement(nodes_coord, natural_coords, shapesd); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template ShapeLagrange::ShapeLagrange(const Mesh & mesh, const ID & id, const MemoryID & memory_id) : ShapeLagrangeBase(mesh, kind, id, memory_id) {} /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_points = integration_points.cols(); UInt nb_element = mesh.getConnectivity(type, ghost_type).size(); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); AKANTU_DEBUG_ASSERT(shape_derivatives.getNbComponent() == size_of_shapesd, "The shapes_derivatives array does not have the correct " << "number of component"); shape_derivatives.resize(nb_element * nb_points); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type, filter_elements); Real * shapesd_val = shape_derivatives.storage(); Array::matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); if (filter_elements != empty_filter) nb_element = filter_elements.size(); for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) { if (filter_elements != empty_filter) shapesd_val = shape_derivatives.storage() + filter_elements(elem) * size_of_shapesd * nb_points; Matrix & X = *x_it; Tensor3 B(shapesd_val, spatial_dimension, nb_nodes_per_element, nb_points); computeShapeDerivativesOnCPointsByElement(X, integration_points, B); if (filter_elements == empty_filter) shapesd_val += size_of_shapesd * nb_points; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange::computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, const ElementType & type, const GhostType & ghost_type, const Array & filter_elements) const { #define AKANTU_COMPUTE_SHAPES(type) \ computeShapeDerivativesOnIntegrationPoints( \ nodes, integration_points, shape_derivatives, ghost_type, \ filter_elements); AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES); #undef AKANTU_COMPUTE_SHAPES } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::precomputeShapesOnIntegrationPoints( const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; Matrix & natural_coords = integration_points(type, ghost_type); UInt size_of_shapes = ElementClass::getShapeSize(); Array & shapes_tmp = shapes.alloc(0, size_of_shapes, itp_type, ghost_type); this->computeShapesOnIntegrationPoints(nodes, natural_coords, shapes_tmp, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::precomputeShapeDerivativesOnIntegrationPoints( const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; Matrix & natural_coords = integration_points(type, ghost_type); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); Array & shapes_derivatives_tmp = shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type); this->computeShapeDerivativesOnIntegrationPoints( nodes, natural_coords, shapes_derivatives_tmp, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_degree_of_freedom, const Array & shapes, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->interpolateElementalFieldOnIntegrationPoints( u_el, out_uq, ghost_type, shapes, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT(shapes.exists(itp_type, ghost_type), "No shapes for the type " << shapes.printType(itp_type, ghost_type)); this->interpolateOnIntegrationPoints(in_u, out_uq, nb_degree_of_freedom, shapes(itp_type, ghost_type), ghost_type, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::gradientOnIntegrationPoints( const Array & in_u, Array & out_nablauq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT( shapes_derivatives.exists(itp_type, ghost_type), "No shapes derivatives for the type " << shapes_derivatives.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->gradientElementalFieldOnIntegrationPoints( u_el, out_nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeBtD( const Array & Ds, Array & BtDs, GhostType ghost_type, const Array & filter_elements) const { auto itp_type = ElementClassProperty::interpolation_type; const auto & shapes_derivatives = this->shapes_derivatives(itp_type, ghost_type); auto spatial_dimension = mesh.getSpatialDimension(); auto size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array shapes_derivatives_filtered(0, shapes_derivatives.getNbComponent()); auto && view = make_view(shapes_derivatives, spatial_dimension, nb_nodes_per_element); auto B_it = view.begin(); auto B_end = view.end(); if (filter_elements != empty_filter) { FEEngine::filterElementalData(this->mesh, shapes_derivatives, shapes_derivatives_filtered, type, ghost_type, filter_elements); auto && view = make_view(shapes_derivatives_filtered, spatial_dimension, nb_nodes_per_element); B_it = view.begin(); B_end = view.end(); } for (auto && values : zip(range(B_it, B_end), - make_view(Ds, spatial_dimension, - Ds.getNbComponent() / spatial_dimension), - make_view(BtDs, BtDs.getNbComponent() / nb_nodes_per_element, nb_nodes_per_element))) { + make_view(Ds, Ds.getNbComponent() / spatial_dimension, + spatial_dimension), + make_view(BtDs, BtDs.getNbComponent() / nb_nodes_per_element, + nb_nodes_per_element))) { const auto & B = std::get<0>(values); const auto & D = std::get<1>(values); auto & Bt_D = std::get<2>(values); // transposed due to the storage layout of B Bt_D.template mul(D, B); } } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeBtDB( const Array & Ds, Array & BtDBs, UInt order_d, GhostType ghost_type, const Array & filter_elements) const { auto itp_type = ElementClassProperty::interpolation_type; const auto & shapes_derivatives = this->shapes_derivatives(itp_type, ghost_type); constexpr auto dim = ElementClass::getSpatialDimension(); auto size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array shapes_derivatives_filtered(0, shapes_derivatives.getNbComponent()); auto && view = make_view(shapes_derivatives, dim, nb_nodes_per_element); auto B_it = view.begin(); auto B_end = view.end(); if (filter_elements != empty_filter) { FEEngine::filterElementalData(this->mesh, shapes_derivatives, shapes_derivatives_filtered, type, ghost_type, filter_elements); auto && view = make_view(shapes_derivatives_filtered, dim, nb_nodes_per_element); B_it = view.begin(); B_end = view.end(); } if (order_d == 4) { UInt tangent_size = VoigtHelper::size; Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); for (auto && values : zip(range(B_it, B_end), make_view(Ds, tangent_size, tangent_size), make_view(BtDBs, dim * nb_nodes_per_element, dim * nb_nodes_per_element))) { const auto & Bfull = std::get<0>(values); const auto & D = std::get<1>(values); auto & Bt_D_B = std::get<2>(values); VoigtHelper::transferBMatrixToSymVoigtBMatrix(Bfull, B, nb_nodes_per_element); Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } } else if (order_d == 2) { Matrix Bt_D(nb_nodes_per_element, dim); - for (auto && values : zip(range(B_it, B_end), make_view(Ds, dim, dim), - make_view(BtDBs, nb_nodes_per_element, - nb_nodes_per_element))) { + for (auto && values : + zip(range(B_it, B_end), make_view(Ds, dim, dim), + make_view(BtDBs, nb_nodes_per_element, nb_nodes_per_element))) { const auto & B = std::get<0>(values); const auto & D = std::get<1>(values); auto & Bt_D_B = std::get<2>(values); Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } } } template <> template <> inline void ShapeLagrange<_ek_regular>::computeBtDB<_point_1>( const Array & /*Ds*/, Array & /*BtDBs*/, UInt /*order_d*/, GhostType /*ghost_type*/, const Array & /*filter_elements*/) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::fieldTimesShapes(const Array & field, Array & field_times_shapes, GhostType ghost_type) const { AKANTU_DEBUG_IN(); field_times_shapes.resize(field.size()); UInt size_of_shapes = ElementClass::getShapeSize(); InterpolationType itp_type = ElementClassProperty::interpolation_type; UInt nb_degree_of_freedom = field.getNbComponent(); const auto & shape = shapes(itp_type, ghost_type); auto field_it = field.begin(nb_degree_of_freedom, 1); auto shapes_it = shape.begin(1, size_of_shapes); auto it = field_times_shapes.begin(nb_degree_of_freedom, size_of_shapes); auto end = field_times_shapes.end(nb_degree_of_freedom, size_of_shapes); for (; it != end; ++it, ++field_it, ++shapes_it) { it->template mul(*field_it, *shapes_it); } AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__ */ diff --git a/src/io/parser/parameter_registry.hh b/src/io/parser/parameter_registry.hh index 56021dff4..76f9e4d7e 100644 --- a/src/io/parser/parameter_registry.hh +++ b/src/io/parser/parameter_registry.hh @@ -1,217 +1,217 @@ /** * @file parameter_registry.hh * * @author Nicolas Richart * * @date creation: Thu Aug 09 2012 * @date last modification: Thu Dec 17 2015 * * @brief Interface of the parameter registry * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "parser.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PARAMETER_REGISTRY_HH__ #define __AKANTU_PARAMETER_REGISTRY_HH__ namespace akantu { class ParserParameter; } namespace akantu { /* -------------------------------------------------------------------------- */ /// Defines the access modes of parsable parameters enum ParameterAccessType { _pat_internal = 0x0001, _pat_writable = 0x0010, _pat_readable = 0x0100, _pat_modifiable = 0x0110, //<_pat_readable | _pat_writable, _pat_parsable = 0x1000, _pat_parsmod = 0x1110 //< _pat_parsable | _pat_modifiable }; /// Bit-wise operator between access modes inline ParameterAccessType operator|(const ParameterAccessType & a, const ParameterAccessType & b) { auto tmp = ParameterAccessType(UInt(a) | UInt(b)); return tmp; } /* -------------------------------------------------------------------------- */ template class ParameterTyped; /** * Interface for the Parameter */ class Parameter { public: Parameter(); Parameter(std::string name, std::string description, ParameterAccessType param_type); virtual ~Parameter() = default; /* ------------------------------------------------------------------------ */ bool isInternal() const; bool isWritable() const; bool isReadable() const; bool isParsable() const; void setAccessType(ParameterAccessType ptype); /* ------------------------------------------------------------------------ */ template void set(const V & value); virtual void setAuto(const ParserParameter & param); template T & get(); template const T & get() const; virtual inline operator Real() const { throw std::bad_cast(); }; template inline operator T() const; /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const; virtual const std::type_info& type() const = 0; protected: /// Returns const instance of templated sub-class ParameterTyped template const ParameterTyped & getParameterTyped() const; /// Returns instance of templated sub-class ParameterTyped template ParameterTyped & getParameterTyped(); protected: /// Name of parameter std::string name; private: /// Description of parameter std::string description; /// Type of access ParameterAccessType param_type{_pat_internal}; }; /* -------------------------------------------------------------------------- */ /* Typed Parameter */ /* -------------------------------------------------------------------------- */ /** * Type parameter transfering a ParserParameter (string: string) to a typed * parameter in the memory of the p */ template class ParameterTyped : public Parameter { public: ParameterTyped(std::string name, std::string description, ParameterAccessType param_type, T & param); /* ------------------------------------------------------------------------ */ template void setTyped(const V & value); void setAuto(const ParserParameter & param) override; T & getTyped(); const T & getTyped() const; void printself(std::ostream & stream) const override; inline operator Real() const override; inline const std::type_info& type() const override { return typeid(T); } private: /// Value of parameter T & param; }; /* -------------------------------------------------------------------------- */ /* Parsable Interface */ /* -------------------------------------------------------------------------- */ /// Defines interface for classes to manipulate parsable parameters class ParameterRegistry { public: ParameterRegistry(); virtual ~ParameterRegistry(); /* ------------------------------------------------------------------------ */ /// Add parameter to the params map template void registerParam(std::string name, T & variable, ParameterAccessType type, const std::string & description = ""); /// Add parameter to the params map (with default value) template void registerParam(std::string name, T & variable, const T & default_value, ParameterAccessType type, const std::string & description = ""); /*------------------------------------------------------------------------- */ protected: void registerSubRegistry(const ID & id, ParameterRegistry & registry); /* ------------------------------------------------------------------------ */ public: /// Set value to a parameter (with possible different type) template void setMixed(const std::string & name, const V & value); /// Set value to a parameter template void set(const std::string & name, const T & value); /// Get value of a parameter inline const Parameter & get(const std::string & name) const; std::vector listParameters() const { std::vector params; for(auto & pair : this->params) params.push_back(pair.first); return params; } std::vector listSubRegisteries() const { std::vector subs; for(auto & pair : this->sub_registries) subs.push_back(pair.first); return subs; } protected: - template T & get(const std::string & name); + template T & get_(const std::string & name); protected: void setParameterAccessType(const std::string & name, ParameterAccessType ptype); /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream, int indent) const; protected: /// Parameters map using Parameters = std::map; Parameters params; /// list of sub-registries using SubRegisteries = std::map; SubRegisteries sub_registries; /// should accessor check in sub registries bool consisder_sub{true}; }; } // akantu #include "parameter_registry_tmpl.hh" #endif /* __AKANTU_PARAMETER_REGISTRY_HH__ */ diff --git a/src/io/parser/parameter_registry_tmpl.hh b/src/io/parser/parameter_registry_tmpl.hh index 3a96f7ad6..935338739 100644 --- a/src/io/parser/parameter_registry_tmpl.hh +++ b/src/io/parser/parameter_registry_tmpl.hh @@ -1,391 +1,391 @@ /** * @file parameter_registry_tmpl.hh * * @author Nicolas Richart * * @date creation: Thu Aug 09 2012 * @date last modification: Fri Mar 27 2015 * * @brief implementation of the templated part of ParameterRegistry class and * the derivated ones * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_error.hh" #include "aka_iterators.hh" #include "parameter_registry.hh" #include "parser.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PARAMETER_REGISTRY_TMPL_HH__ #define __AKANTU_PARAMETER_REGISTRY_TMPL_HH__ namespace akantu { namespace debug { class ParameterException : public Exception { public: ParameterException(const std::string & name, const std::string & message) : Exception(message), name(name) {} const std::string & name; }; class ParameterUnexistingException : public ParameterException { public: ParameterUnexistingException(const std::string & name, const ParameterRegistry & registery) : ParameterException( name, "Parameter " + name + " does not exists in this scope") { auto && params = registery.listParameters(); this->_info = std::accumulate(params.begin(), params.end(), this->_info + "\n Possible parameters are: ", [](auto && str, auto && param) { static auto first = true; auto && ret = str + (first ? " " : ", ") + param; first = false; return ret; }); } }; class ParameterAccessRightException : public ParameterException { public: ParameterAccessRightException(const std::string & name, const std::string & perm) : ParameterException(name, "Parameter " + name + " is not " + perm) {} }; class ParameterWrongTypeException : public ParameterException { public: ParameterWrongTypeException(const std::string & name, const std::type_info & wrong_type, const std::type_info & type) : ParameterException(name, "Parameter " + name + " type error, cannot convert " + debug::demangle(type.name()) + " to " + debug::demangle(wrong_type.name())) {} }; } // namespace debug /* -------------------------------------------------------------------------- */ template const ParameterTyped & Parameter::getParameterTyped() const { try { const auto & tmp = dynamic_cast &>(*this); return tmp; } catch (std::bad_cast &) { AKANTU_CUSTOM_EXCEPTION( debug::ParameterWrongTypeException(name, typeid(T), this->type())); } } /* -------------------------------------------------------------------------- */ template ParameterTyped & Parameter::getParameterTyped() { try { auto & tmp = dynamic_cast &>(*this); return tmp; } catch (std::bad_cast &) { AKANTU_CUSTOM_EXCEPTION( debug::ParameterWrongTypeException(name, typeid(T), this->type())); } } /* ------------------------------------------------------------------------ */ template void Parameter::set(const V & value) { if (!(isWritable())) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "writable")); ParameterTyped & typed_param = getParameterTyped(); typed_param.setTyped(value); } /* ------------------------------------------------------------------------ */ inline void Parameter::setAuto(__attribute__((unused)) const ParserParameter & value) { if (!(isParsable())) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "parsable")); } /* -------------------------------------------------------------------------- */ template const T & Parameter::get() const { if (!(isReadable())) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "readable")); const ParameterTyped & typed_param = getParameterTyped(); return typed_param.getTyped(); } /* -------------------------------------------------------------------------- */ template T & Parameter::get() { ParameterTyped & typed_param = getParameterTyped(); if (!(isReadable()) || !(this->isWritable())) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "accessible")); return typed_param.getTyped(); } /* -------------------------------------------------------------------------- */ template inline Parameter::operator T() const { return this->get(); } /* -------------------------------------------------------------------------- */ template ParameterTyped::ParameterTyped(std::string name, std::string description, ParameterAccessType param_type, T & param) : Parameter(name, description, param_type), param(param) {} /* -------------------------------------------------------------------------- */ template template void ParameterTyped::setTyped(const V & value) { param = value; } /* -------------------------------------------------------------------------- */ template inline void ParameterTyped::setAuto(const ParserParameter & value) { Parameter::setAuto(value); param = static_cast(value); } /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped::setAuto(const ParserParameter & value) { Parameter::setAuto(value); param = value.getValue(); } /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped>::setAuto(const ParserParameter & in_param) { Parameter::setAuto(in_param); Vector tmp = in_param; if (param.size() == 0) { param = tmp; } else { for (UInt i = 0; i < param.size(); ++i) { param(i) = tmp(i); } } } /* -------------------------------------------------------------------------- */ template <> inline void ParameterTyped>::setAuto(const ParserParameter & in_param) { Parameter::setAuto(in_param); Matrix tmp = in_param; if (param.size() == 0) { param = tmp; } else { for (UInt i = 0; i < param.rows(); ++i) { for (UInt j = 0; j < param.cols(); ++j) { param(i, j) = tmp(i, j); } } } } /* -------------------------------------------------------------------------- */ template const T & ParameterTyped::getTyped() const { return param; } /* -------------------------------------------------------------------------- */ template T & ParameterTyped::getTyped() { return param; } /* -------------------------------------------------------------------------- */ template inline void ParameterTyped::printself(std::ostream & stream) const { Parameter::printself(stream); stream << param << "\n"; } /* -------------------------------------------------------------------------- */ template class ParameterTyped> : public Parameter { public: ParameterTyped(std::string name, std::string description, ParameterAccessType param_type, std::vector & param) : Parameter(name, description, param_type), param(param) {} /* ------------------------------------------------------------------------ */ template void setTyped(const V & value) { param = value; } void setAuto(const ParserParameter & value) override { Parameter::setAuto(value); param.clear(); const std::vector & tmp = value; for (auto && z : tmp) { param.emplace_back(z); } } std::vector & getTyped() { return param; } const std::vector & getTyped() const { return param; } void printself(std::ostream & stream) const override { Parameter::printself(stream); stream << "[ "; for (auto && v : param) stream << v << " "; stream << "]\n"; } inline const std::type_info & type() const override { return typeid(std::vector); } private: /// Value of parameter std::vector & param; }; /* ------o-------------------------------------------------------------------- */ template <> inline void ParameterTyped::printself(std::ostream & stream) const { Parameter::printself(stream); stream << std::boolalpha << param << "\n"; } /* -------------------------------------------------------------------------- */ template void ParameterRegistry::registerParam(std::string name, T & variable, ParameterAccessType type, const std::string & description) { auto it = params.find(name); if (it != params.end()) AKANTU_CUSTOM_EXCEPTION(debug::ParameterException( name, "Parameter named " + name + " already registered.")); auto * param = new ParameterTyped(name, description, type, variable); params[name] = param; } /* -------------------------------------------------------------------------- */ template void ParameterRegistry::registerParam(std::string name, T & variable, const T & default_value, ParameterAccessType type, const std::string & description) { variable = default_value; registerParam(name, variable, type, description); } /* -------------------------------------------------------------------------- */ template void ParameterRegistry::setMixed(const std::string & name, const V & value) { auto it = params.find(name); if (it == params.end()) { if (consisder_sub) { for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) { it->second->setMixed(name, value); } } else { AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this)); } } else { Parameter & param = *(it->second); param.set(value); } } /* -------------------------------------------------------------------------- */ template void ParameterRegistry::set(const std::string & name, const T & value) { this->template setMixed(name, value); } /* -------------------------------------------------------------------------- */ -template T & ParameterRegistry::get(const std::string & name) { +template T & ParameterRegistry::get_(const std::string & name) { auto it = params.find(name); if (it == params.end()) { if (consisder_sub) { for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) { try { - return it->second->get(name); + return it->second->get_(name); } catch (...) { } } } // nothing was found not even in sub registries AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this)); } Parameter & param = *(it->second); return param.get(); } /* -------------------------------------------------------------------------- */ const Parameter & ParameterRegistry::get(const std::string & name) const { auto it = params.find(name); if (it == params.end()) { if (consisder_sub) { for (auto it = sub_registries.begin(); it != sub_registries.end(); ++it) { try { return it->second->get(name); } catch (...) { } } } // nothing was found not even in sub registries AKANTU_CUSTOM_EXCEPTION(debug::ParameterUnexistingException(name, *this)); } Parameter & param = *(it->second); return param; } /* -------------------------------------------------------------------------- */ namespace { namespace details { template struct CastHelper { static R convert(const T &) { throw std::bad_cast(); } }; template struct CastHelper::value>> { static R convert(const T & val) { return val; } }; } // namespace details } // namespace template inline ParameterTyped::operator Real() const { if (not isReadable()) AKANTU_CUSTOM_EXCEPTION( debug::ParameterAccessRightException(name, "accessible")); return details::CastHelper::convert(param); } } // namespace akantu #endif /* __AKANTU_PARAMETER_REGISTRY_TMPL_HH__ */ diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh index 0e4b72f0c..c6fc7a798 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,436 +1,437 @@ /** * @file element_type_map.hh * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Fri Oct 02 2015 * * @brief storage class by element type * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_memory.hh" #include "aka_named_argument.hh" #include "element.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_TYPE_MAP_HH__ #define __AKANTU_ELEMENT_TYPE_MAP_HH__ namespace akantu { class FEEngine; } // namespace akantu namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(all_ghost_types); DECLARE_NAMED_ARGUMENT(default_value); DECLARE_NAMED_ARGUMENT(element_kind); DECLARE_NAMED_ARGUMENT(ghost_type); DECLARE_NAMED_ARGUMENT(nb_component); DECLARE_NAMED_ARGUMENT(with_nb_element); DECLARE_NAMED_ARGUMENT(with_nb_nodes_per_element); DECLARE_NAMED_ARGUMENT(spatial_dimension); DECLARE_NAMED_ARGUMENT(do_not_default); } // namespace template class ElementTypeMap; /* -------------------------------------------------------------------------- */ /* ElementTypeMapBase */ /* -------------------------------------------------------------------------- */ /// Common non templated base class for the ElementTypeMap class class ElementTypeMapBase { public: virtual ~ElementTypeMapBase() = default; }; /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template class ElementTypeMap : public ElementTypeMapBase { public: ElementTypeMap(); ~ElementTypeMap() override; inline static std::string printType(const SupportType & type, const GhostType & ghost_type); /*! Tests whether a type is present in the object * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return true if the type is present. */ inline bool exists(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline const Stored & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline Stored & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. THIS METHOD IS * NOT ARRAY SAFE, when using ElementTypeMapArray, use setArray instead * @param data to insert * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline Stored & operator()(const Stored & insert, const SupportType & type, const GhostType & ghost_type = _not_ghost); /// print helper virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ /*! iterator allows to iterate over type-data pairs of the map. The interface * expects the SupportType to be ElementType. */ using DataMap = std::map; class type_iterator : private std::iterator { public: using value_type = const SupportType; using pointer = const SupportType *; using reference = const SupportType &; protected: using DataMapIterator = typename ElementTypeMap::DataMap::const_iterator; public: type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek); type_iterator(const type_iterator & it); type_iterator() = default; inline reference operator*(); inline reference operator*() const; inline type_iterator & operator++(); type_iterator operator++(int); inline bool operator==(const type_iterator & other) const; inline bool operator!=(const type_iterator & other) const; type_iterator & operator=(const type_iterator & other); private: DataMapIterator list_begin; DataMapIterator list_end; UInt dim; ElementKind kind; }; /// helper class to use in range for constructions class ElementTypesIteratorHelper { public: using Container = ElementTypeMap; using iterator = typename Container::type_iterator; ElementTypesIteratorHelper(const Container & container, UInt dim, GhostType ghost_type, ElementKind kind) : container(std::cref(container)), dim(dim), ghost_type(ghost_type), kind(kind) {} template ElementTypesIteratorHelper(const Container & container, use_named_args_t, pack &&... _pack) : ElementTypesIteratorHelper( container, OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), OPTIONAL_NAMED_ARG(ghost_type, _not_ghost), OPTIONAL_NAMED_ARG(element_kind, _ek_regular)) {} ElementTypesIteratorHelper(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(ElementTypesIteratorHelper &&) = default; iterator begin() { return container.get().firstType(dim, ghost_type, kind); } iterator end() { return container.get().lastType(dim, ghost_type, kind); } private: std::reference_wrapper container; UInt dim; GhostType ghost_type; ElementKind kind; }; private: ElementTypesIteratorHelper elementTypesImpl(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const; template ElementTypesIteratorHelper elementTypesImpl(const use_named_args_t & /*unused*/, pack &&... _pack) const; public: template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(use_named_args, std::forward(_pack)...); } template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(std::forward(_pack)...); } public: /*! Get an iterator to the beginning of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the first stored data matching the filters * or an iterator to the end of the map if none match*/ inline type_iterator firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; /*! Get an iterator to the end of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the last stored data matching the filters * or an iterator to the end of the map if none match */ inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; protected: /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline DataMap & getData(GhostType ghost_type); /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline const DataMap & getData(GhostType ghost_type) const; /* ------------------------------------------------------------------------ */ protected: DataMap data; DataMap ghost_data; }; /* -------------------------------------------------------------------------- */ /* Some typedefs */ /* -------------------------------------------------------------------------- */ template class ElementTypeMapArray : public ElementTypeMap *, SupportType>, public Memory { public: using type = T; using array_type = Array; protected: using parent = ElementTypeMap *, SupportType>; using DataMap = typename parent::DataMap; public: using type_iterator = typename parent::type_iterator; /// standard assigment (copy) operator void operator=(const ElementTypeMapArray &) = delete; ElementTypeMapArray(const ElementTypeMapArray &) = delete; /*! Constructor * @param id optional: identifier (string) * @param parent_id optional: parent identifier. for organizational purposes * only * @param memory_id optional: choose a specific memory, defaults to memory 0 */ ElementTypeMapArray(const ID & id = "by_element_type_array", const ID & parent_id = "no_parent", const MemoryID & memory_id = 0) : parent(), Memory(parent_id + ":" + id, memory_id), name(id){}; /*! allocate memory for a new array * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map * @param ghost_type whether to add the field to the data map or the * ghost_data map * @return a reference to the allocated array */ inline Array & alloc(UInt size, UInt nb_component, const SupportType & type, const GhostType & ghost_type, const T & default_value = T()); /*! allocate memory for a new array in both the data and the ghost_data map * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map*/ inline void alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value = T()); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a reference to the array */ inline const Array & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /// access the data of an element, this combine the map and array accessor inline const T & operator()(const Element & element, UInt component = 0) const; /// access the data of an element, this combine the map and array accessor inline T & operator()(const Element & element, UInt component = 0); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a const reference to the array */ inline Array & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @param vect the vector to include into the map * @return stored data corresponding to type. */ inline void setArray(const SupportType & type, const GhostType & ghost_type, const Array & vect); /*! frees all memory related to the data*/ inline void free(); /*! set all values in the ElementTypeMap to zero*/ inline void clear(); /*! set all values in the ElementTypeMap to value */ - inline void set(const T & value); + template + inline void set(const ST & value); /*! deletes and reorders entries in the stored arrays * @param new_numbering a ElementTypeMapArray of new indices. UInt(-1) * indicates * deleted entries. */ inline void onElementsRemoved(const ElementTypeMapArray & new_numbering); /// text output helper void printself(std::ostream & stream, int indent = 0) const override; /*! set the id * @param id the new name */ inline void setID(const ID & id) { this->id = id; } ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const { ElementTypeMap nb_components; for (auto & type : this->elementTypes(dim, ghost_type, kind)) { UInt nb_comp = (*this)(type, ghost_type).getNbComponent(); nb_components(type, ghost_type) = nb_comp; } return nb_components; } /* ------------------------------------------------------------------------ */ /* more evolved allocators */ /* ------------------------------------------------------------------------ */ public: /// initialize the arrays in accordance to a functor template void initialize(const Func & f, const T & default_value, bool do_not_default); /// initialize with sizes and number of components in accordance of a mesh /// content template void initialize(const Mesh & mesh, pack &&... _pack); /// initialize with sizes and number of components in accordance of a fe /// engine content (aka integration points) template void initialize(const FEEngine & fe_engine, pack &&... _pack); /* ------------------------------------------------------------------------ */ /* Accesssors */ /* ------------------------------------------------------------------------ */ public: /// get the name of the internal field AKANTU_GET_MACRO(Name, name, ID); /// name of the elment type map: e.g. connectivity, grad_u ID name; }; /// to store data Array by element type using ElementTypeMapReal = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapInt = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapUInt = ElementTypeMapArray; /// Map of data of type UInt stored in a mesh using UIntDataMap = std::map *>; using ElementTypeMapUIntDataMap = ElementTypeMap; } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_HH__ */ diff --git a/src/mesh/element_type_map_tmpl.hh b/src/mesh/element_type_map_tmpl.hh index c74978ee8..af2a331d6 100644 --- a/src/mesh/element_type_map_tmpl.hh +++ b/src/mesh/element_type_map_tmpl.hh @@ -1,660 +1,661 @@ /** * @file element_type_map_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Fri Oct 02 2015 * * @brief implementation of template functions of the ElementTypeMap and * ElementTypeMapArray classes * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_static_if.hh" #include "element_type_map.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include "element_type_conversion.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ #define __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template inline std::string ElementTypeMap::printType(const SupportType & type, const GhostType & ghost_type) { std::stringstream sstr; sstr << "(" << ghost_type << ":" << type << ")"; return sstr.str(); } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::exists( const SupportType & type, const GhostType & ghost_type) const { return this->getData(ghost_type).find(type) != this->getData(ghost_type).end(); } /* -------------------------------------------------------------------------- */ template inline const Stored & ElementTypeMap:: operator()(const SupportType & type, const GhostType & ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMap::printType(type, ghost_type) << " in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); return it->second; } /* -------------------------------------------------------------------------- */ template inline Stored & ElementTypeMap:: operator()(const SupportType & type, const GhostType & ghost_type) { return this->getData(ghost_type)[type]; } /* -------------------------------------------------------------------------- */ template inline Stored & ElementTypeMap:: operator()(const Stored & insert, const SupportType & type, const GhostType & ghost_type) { auto it = this->getData(ghost_type).find(type); if (it != this->getData(ghost_type).end()) { AKANTU_SILENT_EXCEPTION("Element of type " << ElementTypeMap::printType(type, ghost_type) << " already in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); } else { auto & data = this->getData(ghost_type); const auto & res = data.insert(std::pair(type, insert)); it = res.first; } return it->second; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) { if (ghost_type == _not_ghost) return data; return ghost_data; } /* -------------------------------------------------------------------------- */ template inline const typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) const { if (ghost_type == _not_ghost) return data; return ghost_data; } /* -------------------------------------------------------------------------- */ /// Works only if stored is a pointer to a class with a printself method template void ElementTypeMap::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> [" << std::endl; for (auto gt : ghost_types) { const DataMap & data = getData(gt); for (auto & pair : data) { stream << space << space << ElementTypeMap::printType(pair.first, gt) << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template ElementTypeMap::ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ template ElementTypeMap::~ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ /* ElementTypeMapArray */ /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray::alloc( UInt size, UInt nb_component, const SupportType & type, const GhostType & ghost_type, const T & default_value) { std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; Array * tmp; auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) { std::stringstream sstr; sstr << this->id << ":" << type << ghost_id; tmp = &(Memory::alloc(sstr.str(), size, nb_component, default_value)); std::stringstream sstrg; sstrg << ghost_type; // tmp->setTag(sstrg.str()); this->getData(ghost_type)[type] = tmp; } else { AKANTU_DEBUG_INFO( "The vector " << this->id << this->printType(type, ghost_type) << " already exists, it is resized instead of allocated."); tmp = it->second; it->second->resize(size); } return *tmp; } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value) { this->alloc(size, nb_component, type, _not_ghost, default_value); this->alloc(size, nb_component, type, _ghost, default_value); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::free() { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & pair : data) { dealloc(pair.second->getID()); } data.clear(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::clear() { for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { vect.second->clear(); } } } /* -------------------------------------------------------------------------- */ template -inline void ElementTypeMapArray::set(const T & value) { +template +inline void ElementTypeMapArray::set(const ST & value) { for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { vect.second->set(value); } } } /* -------------------------------------------------------------------------- */ template inline const Array & ElementTypeMapArray:: operator()(const SupportType & type, const GhostType & ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this const ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class(\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray:: operator()(const SupportType & type, const GhostType & ghost_type) { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class (\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::setArray(const SupportType & type, const GhostType & ghost_type, const Array & vect) { auto it = this->getData(ghost_type).find(type); if (AKANTU_DEBUG_TEST(dblWarning) && it != this->getData(ghost_type).end() && it->second != &vect) { AKANTU_DEBUG_WARNING( "The Array " << this->printType(type, ghost_type) << " is already registred, this call can lead to a memory leak."); } this->getData(ghost_type)[type] = &(const_cast &>(vect)); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::onElementsRemoved( const ElementTypeMapArray & new_numbering) { for (auto gt : ghost_types) { for (auto & type : new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) { auto support_type = convertType(type); if (this->exists(support_type, gt)) { const auto & renumbering = new_numbering(type, gt); if (renumbering.size() == 0) continue; auto & vect = this->operator()(support_type, gt); auto nb_component = vect.getNbComponent(); Array tmp(renumbering.size(), nb_component); UInt new_size = 0; for (UInt i = 0; i < vect.size(); ++i) { UInt new_i = renumbering(i); if (new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component, vect.storage() + i * nb_component, nb_component * sizeof(T)); ++new_size; } } tmp.resize(new_size); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ template void ElementTypeMapArray::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; for (UInt g = _not_ghost; g <= _ghost; ++g) { auto gt = (GhostType)g; const DataMap & data = this->getData(gt); typename DataMap::const_iterator it; for (it = data.begin(); it != data.end(); ++it) { stream << space << space << ElementTypeMapArray::printType(it->first, gt) << " [" << std::endl; it->second->printself(stream, indent + 3); stream << space << space << " ]" << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /* SupportType Iterator */ /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek) : list_begin(list_begin), list_end(list_end), dim(dim), kind(ek) {} /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( const type_iterator & it) : list_begin(it.list_begin), list_end(it.list_end), dim(it.dim), kind(it.kind) {} /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator:: operator=(const type_iterator & it) { if (this != &it) { list_begin = it.list_begin; list_end = it.list_end; dim = it.dim; kind = it.kind; } return *this; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() const { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator::operator++() { ++list_begin; while ((list_begin != list_end) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(list_begin->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(list_begin->first))))) ++list_begin; return *this; } /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator ElementTypeMap::type_iterator::operator++(int) { type_iterator tmp(*this); operator++(); return tmp; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator:: operator==(const type_iterator & other) const { return this->list_begin == other.list_begin; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator:: operator!=(const type_iterator & other) const { return this->list_begin != other.list_begin; } /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::ElementTypesIteratorHelper ElementTypeMap::elementTypesImpl(UInt dim, GhostType ghost_type, ElementKind kind) const { return ElementTypesIteratorHelper(*this, dim, ghost_type, kind); } /* -------------------------------------------------------------------------- */ template template typename ElementTypeMap::ElementTypesIteratorHelper ElementTypeMap::elementTypesImpl( const use_named_args_t & unused, pack &&... _pack) const { return ElementTypesIteratorHelper(*this, unused, _pack...); } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator ElementTypeMap::firstType(UInt dim, GhostType ghost_type, ElementKind kind) const { typename DataMap::const_iterator b, e; b = getData(ghost_type).begin(); e = getData(ghost_type).end(); // loop until the first valid type while ((b != e) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(b->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(b->first))))) ++b; return typename ElementTypeMap::type_iterator(b, e, dim, kind); } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator ElementTypeMap::lastType(UInt dim, GhostType ghost_type, ElementKind kind) const { typename DataMap::const_iterator e; e = getData(ghost_type).end(); return typename ElementTypeMap::type_iterator(e, e, dim, kind); } /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const ElementTypeMap & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ class ElementTypeMapArrayInializer { public: ElementTypeMapArrayInializer(UInt spatial_dimension = _all_dimensions, UInt nb_component = 1, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular) : spatial_dimension(spatial_dimension), nb_component(nb_component), ghost_type(ghost_type), element_kind(element_kind) {} const GhostType & ghostType() const { return ghost_type; } protected: UInt spatial_dimension; UInt nb_component; GhostType ghost_type; ElementKind element_kind; }; /* -------------------------------------------------------------------------- */ class MeshElementTypeMapArrayInializer : public ElementTypeMapArrayInializer { public: MeshElementTypeMapArrayInializer( const Mesh & mesh, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular, bool with_nb_element = false, bool with_nb_nodes_per_element = false) : ElementTypeMapArrayInializer(spatial_dimension, nb_component, ghost_type, element_kind), mesh(mesh), with_nb_element(with_nb_element), with_nb_nodes_per_element(with_nb_nodes_per_element) {} decltype(auto) elementTypes() const { return mesh.elementTypes(this->spatial_dimension, this->ghost_type, this->element_kind); } virtual UInt size(const ElementType & type) const { if (with_nb_element) return mesh.getNbElement(type, this->ghost_type); return 0; } UInt nbComponent(const ElementType & type) const { if (with_nb_nodes_per_element) return (this->nb_component * mesh.getNbNodesPerElement(type)); return this->nb_component; } protected: const Mesh & mesh; bool with_nb_element; bool with_nb_nodes_per_element; }; /* -------------------------------------------------------------------------- */ class FEEngineElementTypeMapArrayInializer : public MeshElementTypeMapArrayInializer { public: FEEngineElementTypeMapArrayInializer( const FEEngine & fe_engine, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular); UInt size(const ElementType & type) const override; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; ElementTypesIteratorHelper elementTypes() const; protected: const FEEngine & fe_engine; }; /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const Func & f, const T & default_value, bool do_not_default) { for (auto & type : f.elementTypes()) { if (not this->exists(type, f.ghostType())) if (do_not_default) { auto & array = this->alloc(0, f.nbComponent(type), type, f.ghostType()); array.resize(f.size(type)); } else { this->alloc(f.size(type), f.nbComponent(type), type, f.ghostType(), default_value); } else { auto & array = this->operator()(type, f.ghostType()); if (not do_not_default) array.resize(f.size(type), default_value); else array.resize(f.size(type)); } } } /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const Mesh & mesh, pack &&... _pack) { GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _casper); bool all_ghost_types = requested_ghost_type == _casper; for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; this->initialize( MeshElementTypeMapArrayInializer( mesh, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, mesh.getSpatialDimension()), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_regular), OPTIONAL_NAMED_ARG(with_nb_element, false), OPTIONAL_NAMED_ARG(with_nb_nodes_per_element, false)), OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ template template void ElementTypeMapArray::initialize(const FEEngine & fe_engine, pack &&... _pack) { bool all_ghost_types = OPTIONAL_NAMED_ARG(all_ghost_types, true); GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _not_ghost); for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; this->initialize(FEEngineElementTypeMapArrayInializer( fe_engine, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, UInt(-2)), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_regular)), OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ template inline T & ElementTypeMapArray:: operator()(const Element & element, UInt component) { return this->operator()(element.type, element.ghost_type)(element.element, component); } /* -------------------------------------------------------------------------- */ template inline const T & ElementTypeMapArray:: operator()(const Element & element, UInt component) const { return this->operator()(element.type, element.ghost_type)(element.element, component); } } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ */ diff --git a/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index 2b69e3a2b..f034db4a5 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,783 +1,776 @@ /** * @file heat_transfer_model.cc * * @author Guillaume Anciaux * @author Lucas Frerot * @author David Simon Kammer * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Mon Nov 30 2015 * * @brief Implementation of HeatTransferModel class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "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.cc" #include "mesh.hh" #include "parser.hh" +#include "integrator_gauss.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 &) const { rho.set(model.getCapacity()); } private: const HeatTransferModel & model; }; } } + /* -------------------------------------------------------------------------- */ HeatTransferModel::HeatTransferModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, ModelType::_heat_transfer_model, dim, id, memory_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), conductivity(this->spatial_dimension, this->spatial_dimension) { AKANTU_DEBUG_IN(); this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, _gst_htm_capacity); this->registerSynchronizer(synchronizer, _gst_htm_temperature); this->registerSynchronizer(synchronizer, _gst_htm_gradient_temperature); } - std::stringstream sstr; - sstr << id << ":fem"; - registerFEEngineObject(sstr.str(), mesh, spatial_dimension); + registerFEEngineObject(id + ":fem", mesh, spatial_dimension); #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("paraview_all", 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 dynamic_cast( - getFEEngineClassBoundary(name)); + getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::allocNodalField(Array *& array, const ID & name) { if (array == nullptr) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":" << name; array = &(alloc(sstr_disp.str(), nb_nodes, 1, T())); } } /* -------------------------------------------------------------------------- */ HeatTransferModel::~HeatTransferModel() = default; /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); - auto & fem = getFEEngineClass(); + auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor compute_rho(*this); for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { - fem.assembleFieldLumped(compute_rho, "M", "temperature", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType HeatTransferModel::getMatrixType(const ID & matrix_id) { if (matrix_id == "K" or matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleConductivityMatrix(); } else if (matrix_id == "M") { this->assembleCapacity(); } else { AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { this->assembleCapacityLumped(); } else { AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped() { AKANTU_DEBUG_IN(); if (!this->getDOFManager().hasLumpedMatrix("M")) { this->getDOFManager().getNewLumpedMatrix("M"); } this->getDOFManager().clearLumpedMatrix("M"); assembleCapacityLumped(_not_ghost); assembleCapacityLumped(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->temperature, "temperature"); this->allocNodalField(this->external_heat_rate, "external_heat_rate"); this->allocNodalField(this->internal_heat_rate, "internal_heat_rate"); this->allocNodalField(this->blocked_dofs, "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 == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(this->temperature_rate, "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", _tsst_dynamic_lumped); } case _static: { return std::make_tuple("static", _tsst_static); } case _implicit_dynamic: { return std::make_tuple("implicit", _tsst_dynamic); } default: return std::make_tuple("unknown", _tsst_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions HeatTransferModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic_lumped: { options.non_linear_solver_type = _nls_lumped; options.integration_scheme_type["temperature"] = _ist_forward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["temperature"] = _ist_pseudo_time; options.solution_type["temperature"] = IntegrationScheme::_not_defined; break; } case _tsst_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["temperature"] = _ist_forward_euler; - options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; + options.solution_type["temperature"] = + IntegrationScheme::_temperature_rate; } else { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["temperature"] = _ist_trapezoidal_rule_1; options.solution_type["temperature"] = IntegrationScheme::_temperature; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleConductivityMatrix(bool compute_conductivity) { AKANTU_DEBUG_IN(); if (!this->getDOFManager().hasLumpedMatrix("K")) { this->getDOFManager().getNewLumpedMatrix("K"); } this->getDOFManager().clearLumpedMatrix("K"); switch (mesh.getSpatialDimension()) { case 1: this->assembleConductivityMatrix<1>(_not_ghost, compute_conductivity); break; case 2: this->assembleConductivityMatrix<2>(_not_ghost, compute_conductivity); break; case 3: this->assembleConductivityMatrix<3>(_not_ghost, compute_conductivity); break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::assembleConductivityMatrix(const GhostType & ghost_type, bool compute_conductivity) { AKANTU_DEBUG_IN(); - Mesh & mesh = this->getFEEngine().getMesh(); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type)) { this->assembleConductivityMatrix(type, ghost_type, compute_conductivity); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::assembleConductivityMatrix(const ElementType & type, const GhostType & ghost_type, bool compute_conductivity) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine(); // const Array & shapes_derivatives = // this->getFEEngine().getShapesDerivatives(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array * bt_d_b = new Array(nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); if (compute_conductivity) this->computeConductivityOnQuadPoints(ghost_type); fem.computeBtDB(conductivity_on_qpoints(type, ghost_type), *bt_d_b, 2, type, ghost_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, ghost_type); delete bt_d_b; this->getDOFManager().assembleElementalMatricesToMatrix( "K", "temperature", *K_e, type, ghost_type, _symmetric); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeConductivityOnQuadPoints( const GhostType & ghost_type) { for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { - Array & temperature_interpolated = - temperature_on_qpoints(type, ghost_type); + 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 C_it = conductivity_on_qpoints(type, ghost_type) - .begin(spatial_dimension, spatial_dimension); - auto C_end = conductivity_on_qpoints(type, ghost_type) - .end(spatial_dimension, spatial_dimension); - - auto T_it = temperature_interpolated.begin(); - - for (; C_it != C_end; ++C_it, ++T_it) { - Matrix & C = *C_it; - Real & T = *T_it; + 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; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeKgradT(const GhostType & ghost_type) { - if (this->compute_conductivity) computeConductivityOnQuadPoints(ghost_type); for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { 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->synchronize(_gst_htm_temperature); auto & fem = this->getFEEngine(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(spatial_dimension, ghost_type)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); // compute k \grad T computeKgradT(ghost_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)) { - UInt nb_nodes_per_element = - getFEEngine().getMesh().getNbNodesPerElement(type); + UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array coord(0, nb_nodes_per_element * spatial_dimension); - FEEngine::extractNodalToElementField(getFEEngine().getMesh(), - getFEEngine().getMesh().getNodes(), - coord, type, _not_ghost); + FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, + _not_ghost); auto el_coord = coord.begin(spatial_dimension, nb_nodes_per_element); - UInt nb_element = getFEEngine().getMesh().getNbElement(type); + 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 * density * capacity / conductivitymax; mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ 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; - auto & fem = getFEEngineClass(); + auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor rho_functor(*this); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type)) { fem.assembleFieldMatrix(rho_functor, "M", "temperature", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeRho(Array & rho, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); FEEngine & fem = this->getFEEngine(); - UInt nb_element = fem.getMesh().getNbElement(type, ghost_type); + UInt nb_element = mesh.getNbElement(type, ghost_type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); rho.resize(nb_element * nb_quadrature_points); - Real * rho_1_val = rho.storage(); + rho.set(this->capacity); - /// 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; - } - } + // 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 is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_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, 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 = 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; - Mesh & mesh = getFEEngine().getMesh(); - Mesh::type_iterator it = mesh.firstType(spatial_dimension); - Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); - for (; it != last_type; ++it) { - UInt nb_element = getFEEngine().getMesh().getNbElement(*it, _not_ghost); - UInt nb_quadrature_points = - getFEEngine().getNbIntegrationPoints(*it, _not_ghost); + auto & fem = getFEEngine(); + + for (auto && type : mesh.elementTypes(spatial_dimension)) { + 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 T_it = this->temperature_on_qpoints(*it).begin(); - auto T_end = this->temperature_on_qpoints(*it).end(); + auto T_it = this->temperature_on_qpoints(type).begin(); + auto T_end = this->temperature_on_qpoints(type).end(); getThermalEnergy(Eth_per_quad.begin(), T_it, T_end); - Eth += getFEEngine().integrate(Eth_per_quad, *it); + 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, const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real energy = 0.; if (id == "thermal") energy = getThermalEnergy(type, index); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * HeatTransferModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; dumper::Field * field = nullptr; field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * HeatTransferModel::createNodalFieldReal( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["temperature"] = temperature; real_nodal_fields["temperature_rate"] = temperature_rate; real_nodal_fields["external_heat_rate"] = external_heat_rate; real_nodal_fields["internal_heat_rate"] = internal_heat_rate; real_nodal_fields["capacity_lumped"] = capacity_lumped; real_nodal_fields["increment"] = increment; dumper::Field * field = mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * HeatTransferModel::createElementalField( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const UInt & spatial_dimension, const ElementKind & element_kind) { dumper::Field * field = nullptr; if (field_name == "partitions") field = mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); else if (field_name == "temperature_gradient") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(temperature_gradient, element_kind); field = mesh.createElementalField( temperature_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "conductivity") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(conductivity_on_qpoints, element_kind); field = mesh.createElementalField( conductivity_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ dumper::Field * HeatTransferModel::createElementalField( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const ElementKind & element_kind) { return nullptr; } /* -------------------------------------------------------------------------- */ dumper::Field * HeatTransferModel::createNodalFieldBool( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } /* -------------------------------------------------------------------------- */ dumper::Field * HeatTransferModel::createNodalFieldReal( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } #endif } // akantu diff --git a/src/model/heat_transfer/heat_transfer_model.hh b/src/model/heat_transfer/heat_transfer_model.hh index cb841e85e..29a78b1c2 100644 --- a/src/model/heat_transfer/heat_transfer_model.hh +++ b/src/model/heat_transfer/heat_transfer_model.hh @@ -1,342 +1,341 @@ /** * @file heat_transfer_model.hh * * @author Guillaume Anciaux * @author Lucas Frerot * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Tue Dec 08 2015 * * @brief Model of Heat Transfer * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ - +/* -------------------------------------------------------------------------- */ +#include "data_accessor.hh" +#include "model.hh" +#include "fe_engine.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_HEAT_TRANSFER_MODEL_HH__ #define __AKANTU_HEAT_TRANSFER_MODEL_HH__ -/* -------------------------------------------------------------------------- */ -#include "data_accessor.hh" -#include "integrator_gauss.hh" -#include "model.hh" -#include "shape_lagrange.hh" - namespace akantu { -class IntegrationScheme1stOrder; +template +class IntegratorGauss; +template class ShapeLagrange; } namespace akantu { class HeatTransferModel : public Model, public DataAccessor, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - typedef FEEngineTemplate MyFEEngineType; + using FEEngineType = FEEngineTemplate; HeatTransferModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "heat_transfer_model", const MemoryID & memory_id = 0); virtual ~HeatTransferModel(); /* ------------------------------------------------------------------------ */ /* 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, NonLinearSolverType) override; /// initialize the model void initModel() override; /// allocate all vectors void assembleJacobian(); /// compute the heat flux void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID &) override; /// callback to assemble a Matrix void assembleMatrix(const ID &) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID &) override; std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /* ------------------------------------------------------------------------ */ /* Methods for explicit */ /* ------------------------------------------------------------------------ */ public: /// compute and get the stable time step Real getStableTimeStep(); /// compute the internal heat flux void assembleInternalHeatRate(); /// calculate the lumped capacity vector for heat transfer problem void assembleCapacityLumped(); /* ------------------------------------------------------------------------ */ /* Methods for implicit */ /* ------------------------------------------------------------------------ */ public: /// assemble the conductivity matrix void assembleConductivityMatrix(bool compute_conductivity = true); /// assemble the conductivity matrix void assembleCapacity(); /// assemble the conductivity matrix void assembleCapacity(GhostType ghost_type); /// compute the capacity on quadrature points void computeRho(Array & rho, ElementType type, GhostType ghost_type); protected: /// computation of the residual for the convergence loop void updateResidualInternal(); private: /// compute the heat flux on ghost types void updateResidual(const GhostType & ghost_type, bool compute_conductivity = false); /// calculate the lumped capacity vector for heat transfer problem (w /// ghost type) void assembleCapacityLumped(const GhostType & ghost_type); /// assemble the conductivity matrix (w/ ghost type) template void assembleConductivityMatrix(const GhostType & ghost_type, bool compute_conductivity = true); /// assemble the conductivity matrix template void assembleConductivityMatrix(const ElementType & type, const GhostType & ghost_type, bool compute_conductivity = true); /// compute the conductivity tensor for each quadrature point in an array void computeConductivityOnQuadPoints(const GhostType & ghost_type); /// compute vector k \grad T for each quadrature point void computeKgradT(const GhostType & ghost_type); /// 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: dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Density, density, Real); AKANTU_GET_MACRO(Capacity, capacity, Real); /// get the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// get the assembled heat flux AKANTU_GET_MACRO(InternalHeatRate, *internal_heat_rate, Array &); /// get the lumped capacity AKANTU_GET_MACRO(CapacityLumped, *capacity_lumped, 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(KGradtOnQpoints, k_gradt_on_qpoints, + 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, const ElementType & type, UInt index); /// get the energy denominated by thermal Real getEnergy(const std::string & energy_id); /// get the thermal energy for a given element Real getThermalEnergy(const ElementType & type, UInt index); /// get the thermal energy for a given element Real getThermalEnergy(); protected: /* ------------------------------------------------------------------------ */ FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ----------------------------------------------------------------------- */ template void getThermalEnergy(iterator Eth, Array::const_iterator T_it, Array::const_iterator T_end) const; template void allocNodalField(Array *& array, const ID & name); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// number of iterations UInt n_iter; /// time step Real time_step; /// temperatures array Array * temperature{nullptr}; /// temperatures derivatives array Array * temperature_rate{nullptr}; /// increment array (@f$\delta \dot T@f$ or @f$\delta T@f$) Array * increment{nullptr}; /// the density Real density; /// the speed of the changing temperature ElementTypeMapArray temperature_gradient; /// temperature field on quadrature points ElementTypeMapArray temperature_on_qpoints; /// conductivity tensor on quadrature points ElementTypeMapArray conductivity_on_qpoints; /// vector k \grad T on quad points ElementTypeMapArray k_gradt_on_qpoints; /// external flux vector Array * external_heat_rate{nullptr}; /// residuals array Array * internal_heat_rate{nullptr}; // lumped vector Array * capacity_lumped{nullptr}; /// boundary vector Array * blocked_dofs{nullptr}; // 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; /// analysis method AnalysisMethod method; bool compute_conductivity; }; } // akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_model_inline_impl.cc" #endif /* __AKANTU_HEAT_TRANSFER_MODEL_HH__ */