diff --git a/src/geometry/mesh_geom_common.hh b/src/geometry/mesh_geom_common.hh index 0255a8a74..4d9dbf7cc 100644 --- a/src/geometry/mesh_geom_common.hh +++ b/src/geometry/mesh_geom_common.hh @@ -1,49 +1,59 @@ /** * @file mesh_geom_common.hh * * @author Lucas Frerot * * @date creation: Wed May 13 2015 * @date last modification: Wed May 13 2015 * * @brief Common file for MeshGeom module * * @section LICENSE * * Copyright (©) 2010-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 . * */ #ifndef __AKANTU_MESH_GEOM_COMMON_HH__ #define __AKANTU_MESH_GEOM_COMMON_HH__ #include "aka_common.hh" #include #include #include #include __BEGIN_AKANTU__ +/** + * Cartesian kernel definition + * - CGAL::Cartesian uses internal referene counting + * - CGAL::Simple_cartesian doesn't (better performance) + * Both use exact predicates and inexact constructions (single/double precision) + */ typedef CGAL::Simple_cartesian Cartesian; +/** + * Spherical kernel definition + * Uses exact predicates and inexact constructions + */ typedef CGAL::Spherical_kernel_3 > Spherical; __END_AKANTU__ #endif // __AKANTU_MESH_GEOM_COMMON_HH__ diff --git a/src/geometry/mesh_geom_factory_tmpl.hh b/src/geometry/mesh_geom_factory_tmpl.hh index b7fb524a4..12f1b9aef 100644 --- a/src/geometry/mesh_geom_factory_tmpl.hh +++ b/src/geometry/mesh_geom_factory_tmpl.hh @@ -1,185 +1,187 @@ /** * @file mesh_geom_factory_tmpl.hh * * @author Lucas Frérot * * @date creation: Thu Feb 26 2015 * @date last modification: Fri Mar 6 2015 * * @brief Class for constructing the CGAL primitives of a mesh * * @section LICENSE * * Copyright (©) 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 . * */ #ifndef __AKANTU_MESH_GEOM_FACTORY_TMPL_HH__ #define __AKANTU_MESH_GEOM_FACTORY_TMPL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh_geom_common.hh" #include "mesh_geom_factory.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ template MeshGeomFactory::MeshGeomFactory(const Mesh & mesh) : MeshGeomAbstract(mesh), data_tree(NULL), primitive_list() {} template MeshGeomFactory::~MeshGeomFactory() { delete data_tree; } /** * This function loops over the elements of `type` in the mesh and creates the * AABB tree of geometrical primitves (`data_tree`). */ template void MeshGeomFactory::constructData() { AKANTU_DEBUG_IN(); const GhostType ghost_type = _not_ghost; primitive_list.clear(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type, ghost_type); const Array & nodes = mesh.getNodes(); Array::const_vector_iterator begin = connectivity.begin(nb_nodes_per_element); Array::const_vector_iterator it = connectivity.begin(nb_nodes_per_element); Array::const_vector_iterator end = connectivity.end(nb_nodes_per_element); // This loop builds the list of primitives for (; it != end ; ++it) { const Vector & el_connectivity = *it; Matrix node_coordinates(dim, nb_nodes_per_element); for (UInt i = 0 ; i < nb_nodes_per_element ; i++) for (UInt j = 0 ; j < dim ; j++) node_coordinates(j, i) = nodes(el_connectivity(i), j); addPrimitive(node_coordinates, it - begin); } delete data_tree; + // This condition allows the use of the mesh geom module + // even if types are not compatible with AABB tree algorithm if (TreeTypeHelper::is_valid) data_tree = new typename TreeTypeHelper::tree(primitive_list.begin(), primitive_list.end()); AKANTU_DEBUG_OUT(); } template void MeshGeomFactory::addPrimitive(const Matrix & node_coordinates, UInt id) { this->addPrimitive(node_coordinates, id, this->primitive_list); } -// 2D and _triangle_3 implementation +// (2D, _triangle_3) decomposed into Triangle template<> inline void MeshGeomFactory<2, _triangle_3, Triangle, Cartesian>::addPrimitive( const Matrix & node_coordinates, UInt id, TreeTypeHelper, Cartesian>::container_type & list) { TreeTypeHelper, Cartesian>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), 0.), b(node_coordinates(0, 1), node_coordinates(1, 1), 0.), c(node_coordinates(0, 2), node_coordinates(1, 2), 0.); Triangle t(a, b, c); t.setId(id); list.push_back(t); } -// 2D and _triangle_3 with Line_arc implementation +// (2D, _triangle_3) decomposed into Line_arc template<> inline void MeshGeomFactory<2, _triangle_3, Line_arc, Spherical>::addPrimitive( const Matrix & node_coordinates, UInt id, TreeTypeHelper, Spherical>::container_type & list) { TreeTypeHelper, Spherical>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), 0.), b(node_coordinates(0, 1), node_coordinates(1, 1), 0.), c(node_coordinates(0, 2), node_coordinates(1, 2), 0.); /*std::cout << "elem " << id << " node 1 : x_node=" << node_coordinates(0, 0) << ", x_arc_node=" << a.x() << ", y_node=" << node_coordinates(1, 0) << ", y_arc_node=" << a.y() << std::endl ; std::cout << "elem " << id << " node 2 : x_node=" << node_coordinates(0, 1) << ", x_arc_node=" << b.x() << ", y_node=" << node_coordinates(1, 1) << ", y_arc_node=" << b.y() << std::endl ; std::cout << "elem " << id << " node 2 : x_node=" << node_coordinates(0, 2) << ", x_arc_node=" << c.x() << ", y_node=" << node_coordinates(1, 2) << ", y_arc_node=" << c.y() << std::endl ;*/ CGAL::Line_3 l1(a, b), l2(b, c), l3(c, a); Line_arc s1(l1,a, b), s2(l2, b, c), s3(l3, c, a); s1.setId(id); s1.setSegId(1); s2.setId(id); s2.setSegId(2); s3.setId(id); s3.setSegId(3); list.push_back(s1); list.push_back(s2); list.push_back(s3); } -// 3D and _tetrahedron_4 with triangles implementation +// (3D, _tetrahedron_4) decomposed into Triangle template<> inline void MeshGeomFactory<3, _tetrahedron_4, Triangle, Cartesian>::addPrimitive( const Matrix & node_coordinates, UInt id, TreeTypeHelper, Cartesian>::container_type & list) { TreeTypeHelper, Cartesian>::point_type a(node_coordinates(0, 0), node_coordinates(1, 0), node_coordinates(2, 0)), b(node_coordinates(0, 1), node_coordinates(1, 1), node_coordinates(2, 1)), c(node_coordinates(0, 2), node_coordinates(1, 2), node_coordinates(2, 2)), d(node_coordinates(0, 3), node_coordinates(1, 3), node_coordinates(2, 3)); Triangle t1(a, b, c), t2(b, c, d), t3(c, d, a), t4(d, a, b); t1.setId(id); t2.setId(id); t3.setId(id); t4.setId(id); list.push_back(t1); list.push_back(t2); list.push_back(t3); list.push_back(t4); } __END_AKANTU__ #endif // __AKANTU_MESH_GEOM_FACTORY_TMPL_HH__ diff --git a/src/geometry/tree_type_helper.hh b/src/geometry/tree_type_helper.hh index 60f5d9c71..1f5132501 100644 --- a/src/geometry/tree_type_helper.hh +++ b/src/geometry/tree_type_helper.hh @@ -1,110 +1,111 @@ /** * @file tree_type_helper.hh * * @author Lucas Frérot * * @date creation: Fri Feb 27 2015 * @date last modification: Thu Mar 5 2015 * * @brief Converts element types of a mesh to CGAL primitive types * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TREE_TYPE_HELPER_HH__ #define __AKANTU_TREE_TYPE_HELPER_HH__ #include "aka_common.hh" #include "line_arc.hh" #include "triangle.hh" #include "tetrahedron.hh" #include "aabb_primitive.hh" #include "mesh_geom_common.hh" #include #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ +/// Replacement class for algorithm that can't use the AABB tree types template struct VoidTree { VoidTree(const iterator & begin, const iterator & end) {} }; /// Helper class used to ease the use of CGAL AABB tree algorithm template struct TreeTypeHelper { static const bool is_valid = false; typedef Primitive primitive_type; typedef typename std::list container_type; typedef typename container_type::iterator iterator; typedef typename container_type::const_iterator const_iterator; typedef typename CGAL::Point_3 point_type; typedef VoidTree tree; }; /// Helper class used to ease the use of intersections template struct IntersectionTypeHelper; /** * Macro used to specialize TreeTypeHelper * @param my_primitive associated primitive type * @param my_query query_type * @param my_kernel kernel type */ #define TREE_TYPE_HELPER_MACRO(my_primitive, my_query, my_kernel) \ template<> \ struct TreeTypeHelper, my_kernel> { \ static const bool is_valid = true; \ typedef my_primitive primitive_type; \ typedef my_primitive##_primitive aabb_primitive_type; \ typedef CGAL::Point_3 point_type; \ \ typedef std::list container_type; \ typedef container_type::iterator iterator; \ typedef CGAL::AABB_traits aabb_traits_type; \ typedef CGAL::AABB_tree tree; \ typedef tree::Primitive_id id_type; \ }; \ \ template<> \ struct IntersectionTypeHelper, my_kernel>, my_query> { \ typedef boost::optional< \ TreeTypeHelper, my_kernel>::tree::Intersection_and_primitive_id::Type \ > intersection_type; \ } TREE_TYPE_HELPER_MACRO(Triangle, Cartesian::Segment_3, Cartesian); //TREE_TYPE_HELPER_MACRO(Line_arc, Spherical::Sphere_3, Spherical); #undef TREE_TYPE_HELPER_MACRO __END_AKANTU__ #endif // __AKANTU_TREE_TYPE_HELPER_HH__ diff --git a/src/model/solid_mechanics/embedded_interface_intersector.cc b/src/model/solid_mechanics/embedded_interface_intersector.cc index 2af937fb8..319225270 100644 --- a/src/model/solid_mechanics/embedded_interface_intersector.cc +++ b/src/model/solid_mechanics/embedded_interface_intersector.cc @@ -1,153 +1,161 @@ /** * @file embedded_interface_intersector.cc * * @author Lucas Frerot * * @date creation: Wed Apr 29 2015 * @date last modification: Wed Apr 29 2015 * * @brief Class that loads the interface from mesh and computes intersections * * @section LICENSE * * Copyright (©) 2010-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 "embedded_interface_intersector.hh" #include "mesh_segment_intersector.hh" +/// Helper macro for types in the mesh. Creates an intersector and computes intersection queries #define INTERFACE_INTERSECTOR_CASE(dim, type) do { \ MeshSegmentIntersector intersector(this->mesh, interface_mesh); \ name_to_primitives_it = name_to_primitives_map.begin(); \ for (; name_to_primitives_it != name_to_primitives_end ; ++name_to_primitives_it) { \ intersector.computeIntersectionQueryList( \ name_to_primitives_it->second, \ name_to_primitives_it->first); \ } } while(0) #define INTERFACE_INTERSECTOR_CASE_2D(type) INTERFACE_INTERSECTOR_CASE(2, type) #define INTERFACE_INTERSECTOR_CASE_3D(type) INTERFACE_INTERSECTOR_CASE(3, type) __BEGIN_AKANTU__ EmbeddedInterfaceIntersector::EmbeddedInterfaceIntersector(const Mesh & mesh, const Mesh & primitive_mesh) : MeshGeomAbstract(mesh), interface_mesh(mesh.getSpatialDimension(), "interface_mesh"), primitive_mesh(primitive_mesh) { + // Initiating mesh connectivity and data interface_mesh.addConnectivityType(_segment_2, _not_ghost); interface_mesh.addConnectivityType(_segment_2, _ghost); interface_mesh.registerData("associated_element").alloc(0, 1, _segment_2); interface_mesh.registerData("physical_names").alloc(0, 1, _segment_2); } EmbeddedInterfaceIntersector::~EmbeddedInterfaceIntersector() {} void EmbeddedInterfaceIntersector::constructData() { AKANTU_DEBUG_IN(); const UInt dim = this->mesh.getSpatialDimension(); + if (dim == 1) + AKANTU_DEBUG_ERROR("No embedded model in 1D. Deactivate intersection initialization"); + Array * physical_names = NULL; try { physical_names = &const_cast &>(this->primitive_mesh.getData("physical_names", _segment_2)); } catch (debug::Exception & e) { AKANTU_DEBUG_ERROR("You must define physical names to reinforcements in order to use the embedded model"); throw e; } const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_segment_2); Array::const_vector_iterator connectivity = primitive_mesh.getConnectivity(_segment_2).begin(nb_nodes_per_element); Array::scalar_iterator names_it = physical_names->begin(), names_end = physical_names->end(); std::map > name_to_primitives_map; + // Loop over the physical names and register segment lists in name_to_primitives_map for (; names_it != names_end ; ++names_it) { UInt element_id = names_it - physical_names->begin(); const Vector el_connectivity = connectivity[element_id]; K::Segment_3 segment = this->createSegment(el_connectivity); name_to_primitives_map[*names_it].push_back(segment); } // Loop over the background types of the mesh Mesh::type_iterator type_it = this->mesh.firstType(dim, _not_ghost), type_end = this->mesh.lastType(dim, _not_ghost); std::map >::iterator name_to_primitives_it, name_to_primitives_end = name_to_primitives_map.end(); for (; type_it != type_end ; ++type_it) { + // Used in AKANTU_BOOST_ELEMENT_SWITCH ElementType type = *type_it; AKANTU_DEBUG_INFO("Computing intersections with background element type " << type); switch(dim) { case 1: - AKANTU_DEBUG_ERROR("No embedded model in 1D"); break; case 2: + // Compute intersections for supported 2D elements AKANTU_BOOST_ELEMENT_SWITCH(INTERFACE_INTERSECTOR_CASE_2D, (_triangle_3)); break; case 3: + // Compute intersections for supported 3D elements AKANTU_BOOST_ELEMENT_SWITCH(INTERFACE_INTERSECTOR_CASE_3D, (_tetrahedron_4)); break; } } AKANTU_DEBUG_OUT(); } K::Segment_3 EmbeddedInterfaceIntersector::createSegment(const Vector & connectivity) { AKANTU_DEBUG_IN(); K::Point_3 * source = NULL, * target = NULL; const Array & nodes = this->primitive_mesh.getNodes(); if (this->mesh.getSpatialDimension() == 2) { source = new K::Point_3(nodes(connectivity(0), 0), nodes(connectivity(0), 1), 0.); target = new K::Point_3(nodes(connectivity(1), 0), nodes(connectivity(1), 1), 0.); } else if (this->mesh.getSpatialDimension() == 3) { source = new K::Point_3(nodes(connectivity(0), 0), nodes(connectivity(0), 1), nodes(connectivity(0), 2)); target = new K::Point_3(nodes(connectivity(1), 0), nodes(connectivity(1), 1), nodes(connectivity(1), 2)); } K::Segment_3 segment(*source, *target); delete source; delete target; AKANTU_DEBUG_OUT(); return segment; } __END_AKANTU__ #undef INTERFACE_INTERSECTOR_CASE #undef INTERFACE_INTERSECTOR_CASE_2D #undef INTERFACE_INTERSECTOR_CASE_3D diff --git a/src/model/solid_mechanics/embedded_interface_model.cc b/src/model/solid_mechanics/embedded_interface_model.cc index 7b38c6f6b..85b5fd5d1 100644 --- a/src/model/solid_mechanics/embedded_interface_model.cc +++ b/src/model/solid_mechanics/embedded_interface_model.cc @@ -1,157 +1,164 @@ /** * @file embedded_interface_model.cc * * @author Lucas Frérot * * @date creation: Mon Mar 9 2015 * @date last modification: Mon Mar 9 2015 * * @brief Model of Solid Mechanics with embedded interfaces * * @section LICENSE * * Copyright (©) 2010-2012, 2014 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 "embedded_interface_model.hh" #include "material_reinforcement.hh" #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" # include "dumpable_inline_impl.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ const EmbeddedInterfaceModelOptions default_embedded_interface_model_options(_explicit_lumped_mass, false, false); EmbeddedInterfaceModel::EmbeddedInterfaceModel(Mesh & mesh, Mesh & primitive_mesh, UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, spatial_dimension, id, memory_id), intersector(mesh, primitive_mesh), interface_mesh(NULL), primitive_mesh(primitive_mesh), interface_material_selector(NULL) { // This pointer should be deleted by ~SolidMechanicsModel() MaterialSelector * mat_sel_pointer = new MeshDataMaterialSelector("physical_names", *this); this->setMaterialSelector(*mat_sel_pointer); interface_mesh = &(intersector.getInterfaceMesh()); + + // Create 1D FEEngine on the interface mesh registerFEEngineObject("EmbeddedInterfaceFEEngine", *interface_mesh, 1); } EmbeddedInterfaceModel::~EmbeddedInterfaceModel() { delete interface_material_selector; } void EmbeddedInterfaceModel::initFull(const ModelOptions & options) { const EmbeddedInterfaceModelOptions & eim_options = dynamic_cast(options); // We don't want to initiate materials before shape functions are initialized SolidMechanicsModelOptions dummy_options(eim_options.analysis_method, true); + // Do no initialize interface_mesh if told so if (!eim_options.no_init_intersections) intersector.constructData(); + // Initialize interface FEEngine FEEngine & engine = getFEEngine("EmbeddedInterfaceFEEngine"); engine.initShapeFunctions(_not_ghost); engine.initShapeFunctions(_ghost); SolidMechanicsModel::initFull(dummy_options); + // This will call SolidMechanicsModel::initMaterials() last this->initMaterials(); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("reinforcement", id); this->mesh.addDumpMeshToDumper("reinforcement", *interface_mesh, 1, _not_ghost, _ek_regular); #endif } +// This function is very similar to SolidMechanicsModel's void EmbeddedInterfaceModel::initMaterials() { Element element; delete interface_material_selector; interface_material_selector = new InterfaceMeshDataMaterialSelector("physical_names", *this); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { element.ghost_type = *gt; Mesh::type_iterator it = interface_mesh->firstType(1, *gt); Mesh::type_iterator end = interface_mesh->lastType(1, *gt); for (; it != end ; ++it) { UInt nb_element = interface_mesh->getNbElement(*it, *gt); element.type = *it; Array & mat_indexes = material_index.alloc(nb_element, 1, *it, *gt); for (UInt el = 0 ; el < nb_element ; el++) { element.element = el; UInt mat_index = (*interface_material_selector)(element); AKANTU_DEBUG_ASSERT(mat_index < materials.size(), "The material selector returned an index that does not exist"); mat_indexes(element.element) = mat_index; materials.at(mat_index)->addElement(*it, el, *gt); } } } SolidMechanicsModel::initMaterials(); } void EmbeddedInterfaceModel::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { #ifdef AKANTU_USE_IOHELPER dumper::Field * field = NULL; + // If dumper is reinforcement, create a 1D elemental field if (dumper_name == "reinforcement") field = this->createElementalField(field_id, group_name, padding_flag, 1, element_kind); else SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, element_kind, padding_flag); if (field) { DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name,group_name); Model::addDumpGroupFieldToDumper(field_id,field,dumper); } #endif } __END_AKANTU__ diff --git a/src/model/solid_mechanics/embedded_interface_model.hh b/src/model/solid_mechanics/embedded_interface_model.hh index bfaed5e8f..20d6c28f9 100644 --- a/src/model/solid_mechanics/embedded_interface_model.hh +++ b/src/model/solid_mechanics/embedded_interface_model.hh @@ -1,153 +1,163 @@ /** * @file embedded_interface_model.hh * * @author Lucas Frérot * * @date creation: Mon Mar 9 2015 * @date last modification: Mon Mar 9 2015 * * @brief Model of Solid Mechanics with embedded interfaces * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_EMBEDDED_INTERFACE_MODEL_HH__ #define __AKANTU_EMBEDDED_INTERFACE_MODEL_HH__ #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "mesh.hh" #include "embedded_interface_intersector.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /// Options for the EmbeddedInterfaceModel struct EmbeddedInterfaceModelOptions : SolidMechanicsModelOptions { + /** + * @brief Constructor for EmbeddedInterfaceModelOptions + * @param analysis_method see SolidMeechanicsModelOptions + * @param no_init_intersections ignores the embedded elements + * @param no_init_materials see SolidMeechanicsModelOptions + */ EmbeddedInterfaceModelOptions(AnalysisMethod analysis_method = _explicit_lumped_mass, bool no_init_intersections = false, bool no_init_materials = false): SolidMechanicsModelOptions(analysis_method, no_init_materials), no_init_intersections(no_init_intersections) {} + /// Ignore the reinforcements bool no_init_intersections; }; extern const EmbeddedInterfaceModelOptions default_embedded_interface_model_options; /** * @brief Solid mechanics model using the embedded model. * * This SolidMechanicsModel subclass implements the embedded model, * a method used to represent 1D elements in a finite elements model * (eg. reinforcements in concrete). * * In addition to the SolidMechanicsModel properties, this model has * a mesh of the 1D elements embedded in the model, and an instance of the * EmbeddedInterfaceIntersector class for the computation of the intersections of the * 1D elements with the background (bulk) mesh. * * @see MaterialReinforcement */ class EmbeddedInterfaceModel : public SolidMechanicsModel { + /// Same type as SolidMechanicsModel typedef FEEngineTemplate MyFEEngineType; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /** * @brief Constructor * * @param mesh main mesh (concrete) * @param primitive_mesh mesh of the embedded reinforcement */ EmbeddedInterfaceModel(Mesh & mesh, Mesh & primitive_mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "embedded_interface_model", const MemoryID & memory_id = 0); /// Destructor virtual ~EmbeddedInterfaceModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Initialise the model virtual void initFull(const ModelOptions & options = default_embedded_interface_model_options); /// Initialise the materials virtual void initMaterials(); + /// Allows filtering of dump fields which need to be dumpes on interface mesh virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - /// get interface mesh + /// Get interface mesh AKANTU_GET_MACRO(InterfaceMesh, *interface_mesh, Mesh &); - /// get associated elements + /// Get associated elements AKANTU_GET_MACRO_BY_ELEMENT_TYPE(InterfaceAssociatedElements, interface_mesh->getData("associated_element"), Element); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Intersector object to build the interface mesh EmbeddedInterfaceIntersector intersector; /// Interface mesh (weak reference) Mesh * interface_mesh; /// Mesh used to create the CGAL primitives for intersections Mesh & primitive_mesh; /// Material selector for interface MaterialSelector * interface_material_selector; }; +/// Material selector based on mesh data for interface elements template class InterfaceMeshDataMaterialSelector : public ElementDataMaterialSelector { public: InterfaceMeshDataMaterialSelector(const std::string & name, const EmbeddedInterfaceModel & model, UInt first_index = 1) : ElementDataMaterialSelector(model.getInterfaceMesh().getData(name), model, first_index) {} }; __END_AKANTU__ #endif // __AKANTU_EMBEDDED_INTERFACE_MODEL_HH__ diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc index b3b2dde16..a2a008579 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.cc @@ -1,758 +1,752 @@ /** * @file material_reinforcement.cc * * @author Lucas Frérot * * @date creation: Thu Mar 12 2015 * @date last modification: Thu Mar 12 2015 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_voigthelper.hh" #include "material_reinforcement.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialReinforcement::MaterialReinforcement(SolidMechanicsModel & model, UInt spatial_dimension, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Material(model, dim, mesh, fe_engine, id), // /!\ dim, not spatial_dimension ! model(NULL), stress_embedded("stress_embedded", *this, 1, fe_engine, this->element_filter), gradu_embedded("gradu_embedded", *this, 1, fe_engine, this->element_filter), directing_cosines("directing_cosines", *this, 1, fe_engine, this->element_filter), pre_stress("pre_stress", *this, 1, fe_engine, this->element_filter), area(1.0), shape_derivatives() { AKANTU_DEBUG_IN(); this->initialize(model); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initialize(SolidMechanicsModel & a_model) { this->model = dynamic_cast(&a_model); AKANTU_DEBUG_ASSERT(this->model != NULL, "MaterialReinforcement needs an EmbeddedInterfaceModel"); this->registerParam("area", area, _pat_parsable | _pat_modifiable, "Reinforcement cross-sectional area"); this->registerParam("pre_stress", pre_stress, _pat_parsable | _pat_modifiable, "Uniform pre-stress"); + // Reallocate the element filter this->element_filter.free(); this->model->getInterfaceMesh().initElementTypeMapArray(this->element_filter, 1, 1, false, _ek_regular); } /* -------------------------------------------------------------------------- */ template MaterialReinforcement::~MaterialReinforcement() { AKANTU_DEBUG_IN(); ElementTypeMap *>::type_iterator it = shape_derivatives.firstType(); ElementTypeMap *>::type_iterator end = shape_derivatives.lastType(); for (; it != end ; ++it) { delete shape_derivatives(*it, _not_ghost); delete shape_derivatives(*it, _ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initMaterial() { Material::initMaterial(); - stress_embedded.initialize(dim * dim); // Check this + stress_embedded.initialize(dim * dim); gradu_embedded.initialize(dim * dim); pre_stress.initialize(1); /// We initialise the stuff that is not going to change during the simulation this->allocBackgroundShapeDerivatives(); this->initBackgroundShapeDerivatives(); this->initDirectingCosines(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::allocBackgroundShapeDerivatives() { AKANTU_DEBUG_IN(); Mesh & interface_mesh = model->getInterfaceMesh(); Mesh & mesh = model->getMesh(); - ghost_type_t::iterator int_ghost_it = ghost_type_t::begin(); // Loop over interface ghosts for (; int_ghost_it != ghost_type_t::end() ; ++int_ghost_it) { Mesh::type_iterator interface_type_it = interface_mesh.firstType(1, *int_ghost_it); Mesh::type_iterator interface_type_end = interface_mesh.lastType(1, *int_ghost_it); + // Loop over interface types for (; interface_type_it != interface_type_end ; ++interface_type_it) { Mesh::type_iterator background_type_it = mesh.firstType(dim, *int_ghost_it); Mesh::type_iterator background_type_end = mesh.lastType(dim, *int_ghost_it); + // Loop over background types for (; background_type_it != background_type_end ; ++background_type_it) { const ElementType & int_type = *interface_type_it; const ElementType & back_type = *background_type_it; const GhostType & int_ghost = *int_ghost_it; std::string shaped_id = "embedded_shape_derivatives"; if (int_ghost == _ghost) shaped_id += ":ghost"; ElementTypeMapArray * shaped_etma = new ElementTypeMapArray(shaped_id, this->name); UInt nb_points = Mesh::getNbNodesPerElement(back_type); UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(int_type); UInt nb_elements = element_filter(int_type, int_ghost).getSize(); + // Alloc the background ElementTypeMapArray shaped_etma->alloc(nb_elements * nb_quad_points, dim * nb_points, back_type); + // Insert the background ElementTypeMapArray in the interface ElementTypeMap shape_derivatives(shaped_etma, int_type, int_ghost); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initBackgroundShapeDerivatives() { AKANTU_DEBUG_IN(); Mesh & mesh = model->getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, _not_ghost); Mesh::type_iterator type_end = mesh.lastType(dim, _not_ghost); for (; type_it != type_end ; ++type_it) { computeBackgroundShapeDerivatives(*type_it, _not_ghost); //computeBackgroundShapeDerivatives(*type_it, _ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initDirectingCosines() { AKANTU_DEBUG_IN(); Mesh & mesh = model->getInterfaceMesh(); Mesh::type_iterator type_it = mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = mesh.lastType(1, _not_ghost); const UInt voigt_size = getTangentStiffnessVoigtSize(dim); directing_cosines.initialize(voigt_size * voigt_size); for (; type_it != type_end ; ++type_it) { computeDirectingCosines(*type_it, _not_ghost); - computeDirectingCosines(*type_it, _ghost); + //computeDirectingCosines(*type_it, _ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = model->getInterfaceMesh(); Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); for (; type_it != type_end ; ++type_it) { assembleStiffnessMatrix(*type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); computeAllStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = model->getInterfaceMesh(); Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); for (; type_it != type_end ; ++type_it) { assembleResidual(*type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeGradU(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(type); Array & gradu_vec = gradu_embedded(type, ghost_type); Mesh::type_iterator back_it = model->getMesh().firstType(dim, ghost_type); Mesh::type_iterator back_end = model->getMesh().lastType(dim, ghost_type); for (; back_it != back_end ; ++back_it) { UInt nodes_per_background_e = Mesh::getNbNodesPerElement(*back_it); Array & shapesd = shape_derivatives(type, ghost_type)->operator()(*back_it, ghost_type); Array * background_filter = new Array(nb_element, 1, "background_filter"); - filterInterfaceBackgroundElements(*background_filter, *back_it, type, ghost_type, ghost_type); + filterInterfaceBackgroundElements(*background_filter, *back_it, type, ghost_type); Array * disp_per_element = new Array(0, dim * nodes_per_background_e, "disp_elem"); FEEngine::extractNodalToElementField(model->getMesh(), model->getDisplacement(), *disp_per_element, *back_it, ghost_type, *background_filter); Array::matrix_iterator disp_it = disp_per_element->begin(dim, nodes_per_background_e); Array::matrix_iterator disp_end = disp_per_element->end(dim, nodes_per_background_e); Array::matrix_iterator shapes_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator grad_u_it = gradu_vec.begin(dim, dim); for (; disp_it != disp_end ; ++disp_it) { for (UInt i = 0; i < nb_quad_points; i++, ++shapes_it, ++grad_u_it) { Matrix & B = *shapes_it; Matrix & du = *grad_u_it; Matrix & u = *disp_it; du.mul(u, B); } } delete background_filter; delete disp_per_element; } AKANTU_DEBUG_OUT(); } template void MaterialReinforcement::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh::type_iterator it = model->getInterfaceMesh().firstType(); Mesh::type_iterator last_type = model->getInterfaceMesh().lastType(); for(; it != last_type; ++it) { computeGradU(*it, ghost_type); computeStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleResidual(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & mesh = model->getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); for (; type_it != type_end ; ++type_it) { - assembleResidual(type, *type_it, ghost_type, _not_ghost); - //assembleResidual(type, *type_it, ghost_type, _ghost); + assembleResidual(type, *type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleResidual(const ElementType & interface_type, const ElementType & background_type, - GhostType interface_ghost, - GhostType background_ghost) { + GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt voigt_size = getTangentStiffnessVoigtSize(dim); Array & residual = const_cast &>(model->getResidual()); FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); FEEngine & background_engine = model->getFEEngine(); - Array & elem_filter = element_filter(interface_type, interface_ghost); + Array & elem_filter = element_filter(interface_type, ghost_type); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); - UInt nb_quadrature_points = interface_engine.getNbQuadraturePoints(interface_type, interface_ghost); + UInt nb_quadrature_points = interface_engine.getNbQuadraturePoints(interface_type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt back_dof = dim * nodes_per_background_e; - Array & shapesd = shape_derivatives(interface_type, interface_ghost)->operator()(background_type, background_ghost); + Array & shapesd = shape_derivatives(interface_type, ghost_type)->operator()(background_type, ghost_type); Array * integrant = new Array(nb_quadrature_points * nb_element, back_dof, "integrant"); Array::vector_iterator integrant_it = integrant->begin(back_dof); Array::vector_iterator integrant_end = integrant->end(back_dof); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator C_it = - directing_cosines(interface_type, interface_ghost).begin(voigt_size, voigt_size); + directing_cosines(interface_type, ghost_type).begin(voigt_size, voigt_size); Array::matrix_iterator sigma_it = - stress_embedded(interface_type, interface_ghost).begin(dim, dim); + stress_embedded(interface_type, ghost_type).begin(dim, dim); Vector sigma(voigt_size); Matrix Bvoigt(voigt_size, back_dof); Vector Ct_sigma(voigt_size); for (; integrant_it != integrant_end ; ++integrant_it, ++B_it, ++C_it, ++sigma_it) { VoigtHelper::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt, nodes_per_background_e); Matrix & C = *C_it; Vector & BtCt_sigma = *integrant_it; stressTensorToVoigtVector(*sigma_it, sigma); Ct_sigma.mul(C, sigma); BtCt_sigma.mul(Bvoigt, Ct_sigma); BtCt_sigma *= area; } Array * residual_interface = new Array(nb_element, back_dof, "residual_interface"); interface_engine.integrate(*integrant, *residual_interface, back_dof, - interface_type, interface_ghost, + interface_type, ghost_type, elem_filter); delete integrant; Array * background_filter = new Array(nb_element, 1, "background_filter"); filterInterfaceBackgroundElements(*background_filter, background_type, interface_type, - background_ghost, interface_ghost); + ghost_type); background_engine.assembleArray(*residual_interface, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), - dim, background_type, background_ghost, *background_filter, -1.0); + dim, background_type, ghost_type, *background_filter, -1.0); delete residual_interface; delete background_filter; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::filterInterfaceBackgroundElements(Array & filter, const ElementType & type, const ElementType & interface_type, - GhostType ghost_type, - GhostType interface_ghost_type) { + GhostType ghost_type) { AKANTU_DEBUG_IN(); filter.resize(0); filter.clear(); - Array & elements = model->getInterfaceAssociatedElements(interface_type, interface_ghost_type); - Array & elem_filter = element_filter(interface_type, interface_ghost_type); + Array & elements = model->getInterfaceAssociatedElements(interface_type, ghost_type); + Array & elem_filter = element_filter(interface_type, ghost_type); Array::scalar_iterator filter_it = elem_filter.begin(), filter_end = elem_filter.end(); for (; filter_it != filter_end ; ++filter_it) { Element & elem = elements(*filter_it); if (elem.type == type) filter.push_back(elem.element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeDirectingCosines(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = this->model->getInterfaceMesh(); const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const UInt steel_dof = dim * nb_nodes_per_element; const UInt voigt_size = getTangentStiffnessVoigtSize(dim); const UInt nb_quad_points = model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(type, ghost_type); Array node_coordinates(this->element_filter(type, ghost_type).getSize(), steel_dof); this->model->getFEEngine().template extractNodalToElementField(interface_mesh, interface_mesh.getNodes(), node_coordinates, type, ghost_type, this->element_filter(type, ghost_type)); Array::matrix_iterator directing_cosines_it = directing_cosines(type, ghost_type).begin(voigt_size, voigt_size); Array::matrix_iterator node_coordinates_it = node_coordinates.begin(dim, nb_nodes_per_element); Array::matrix_iterator node_coordinates_end = node_coordinates.end(dim, nb_nodes_per_element); for (; node_coordinates_it != node_coordinates_end ; ++node_coordinates_it) { for (UInt i = 0 ; i < nb_quad_points ; i++, ++directing_cosines_it) { Matrix & nodes = *node_coordinates_it; Matrix & cosines = *directing_cosines_it; computeDirectingCosinesOnQuad(nodes, cosines); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & mesh = model->getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); for (; type_it != type_end ; ++type_it) { - assembleStiffnessMatrix(type, *type_it, ghost_type, _not_ghost); - //assembleStiffnessMatrix(type, *type_it, ghost_type, _ghost); + assembleStiffnessMatrix(type, *type_it, ghost_type); } AKANTU_DEBUG_OUT(); } template void MaterialReinforcement::assembleStiffnessMatrix(const ElementType & interface_type, const ElementType & background_type, - GhostType interface_ghost, - GhostType background_ghost) { + GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt voigt_size = getTangentStiffnessVoigtSize(dim); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); FEEngine & background_engine = model->getFEEngine(); FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); - Array & elem_filter = element_filter(interface_type, interface_ghost); - Array & grad_u = gradu_embedded(interface_type, interface_ghost); + Array & elem_filter = element_filter(interface_type, ghost_type); + Array & grad_u = gradu_embedded(interface_type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); - UInt nb_quadrature_points = interface_engine.getNbQuadraturePoints(interface_type, interface_ghost); + UInt nb_quadrature_points = interface_engine.getNbQuadraturePoints(interface_type, ghost_type); UInt back_dof = dim * nodes_per_background_e; UInt integrant_size = back_dof; grad_u.resize(nb_quadrature_points * nb_element); - //need function model->getInterfaceDisplacement() - /*interface_engine.gradientOnQuadraturePoints(model->getInterfaceDisplacement(), - grad_u, dim, interface_type, interface_ghost, - elem_filter);*/ - Array * tangent_moduli = new Array(nb_element * nb_quadrature_points, 1, "interface_tangent_moduli"); tangent_moduli->clear(); - computeTangentModuli(interface_type, *tangent_moduli, interface_ghost); + computeTangentModuli(interface_type, *tangent_moduli, ghost_type); - Array & shapesd = shape_derivatives(interface_type, interface_ghost)->operator()(background_type, background_ghost); + Array & shapesd = shape_derivatives(interface_type, ghost_type)->operator()(background_type, ghost_type); Array * integrant = new Array(nb_element * nb_quadrature_points, integrant_size * integrant_size, "B^t*C^t*D*C*B"); integrant->clear(); /// Temporary matrices for integrant product Matrix Bvoigt(voigt_size, back_dof); Matrix DC(voigt_size, voigt_size); Matrix DCB(voigt_size, back_dof); Matrix CtDCB(voigt_size, back_dof); Array::scalar_iterator D_it = tangent_moduli->begin(); Array::scalar_iterator D_end = tangent_moduli->end(); Array::matrix_iterator C_it = - directing_cosines(interface_type, interface_ghost).begin(voigt_size, voigt_size); + directing_cosines(interface_type, ghost_type).begin(voigt_size, voigt_size); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator integrant_it = integrant->begin(integrant_size, integrant_size); for (; D_it != D_end ; ++D_it, ++C_it, ++B_it, ++integrant_it) { Real & D = *D_it; Matrix & C = *C_it; Matrix & B = *B_it; Matrix & BtCtDCB = *integrant_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix(B, Bvoigt, nodes_per_background_e); DC.clear(); DC(0, 0) = D * area; DC *= C; DCB.mul(DC, Bvoigt); CtDCB.mul(C, DCB); BtCtDCB.mul(Bvoigt, CtDCB); } delete tangent_moduli; Array * K_interface = new Array(nb_element, integrant_size * integrant_size, "K_interface"); interface_engine.integrate(*integrant, *K_interface, integrant_size * integrant_size, - interface_type, interface_ghost, + interface_type, ghost_type, elem_filter); delete integrant; Array * background_filter = new Array(nb_element, 1, "background_filter"); filterInterfaceBackgroundElements(*background_filter, background_type, interface_type, - background_ghost, interface_ghost); + ghost_type); - background_engine.assembleMatrix(*K_interface, K, dim, background_type, background_ghost, *background_filter); + background_engine.assembleMatrix(*K_interface, K, dim, background_type, ghost_type, *background_filter); delete K_interface; delete background_filter; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -/// In this function, type and ghost type refer to background elements +/// In this function, type and ghost_type refer to background elements template void MaterialReinforcement::computeBackgroundShapeDerivatives(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = model->getInterfaceMesh(); FEEngine & engine = model->getFEEngine(); FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); Mesh::type_iterator interface_type = interface_mesh.firstType(); Mesh::type_iterator interface_last = interface_mesh.lastType(); for (; interface_type != interface_last ; ++interface_type) { Array & filter = element_filter(*interface_type, ghost_type); const UInt nb_elements = filter.getSize(); const UInt nb_nodes = Mesh::getNbNodesPerElement(type); const UInt nb_quad_per_element = interface_engine.getNbQuadraturePoints(*interface_type); Array quad_pos(nb_quad_per_element * nb_elements, dim, "interface_quad_points"); quad_pos.resize(nb_quad_per_element * nb_elements); interface_engine.interpolateOnQuadraturePoints(interface_mesh.getNodes(), quad_pos, dim, *interface_type, ghost_type, filter); Array & background_shapesd = shape_derivatives(*interface_type, ghost_type)->operator()(type, ghost_type); background_shapesd.clear(); Array * background_elements = new Array(nb_elements, 1, "computeBackgroundShapeDerivatives:background_filter"); - filterInterfaceBackgroundElements(*background_elements, type, *interface_type, ghost_type, ghost_type); + filterInterfaceBackgroundElements(*background_elements, type, *interface_type, ghost_type); Array::scalar_iterator back_it = background_elements->begin(), back_end = background_elements->end(); Array::matrix_iterator shapesd_it = background_shapesd.begin(dim, nb_nodes); Array::vector_iterator quad_pos_it = quad_pos.begin(dim); for (; back_it != back_end ; ++back_it) { for (UInt i = 0 ; i < nb_quad_per_element ; i++, ++shapesd_it, ++quad_pos_it) engine.computeShapeDerivatives(*quad_pos_it, *back_it, type, *shapesd_it, ghost_type); } delete background_elements; } AKANTU_DEBUG_OUT(); } template Real MaterialReinforcement::getEnergy(std::string id) { AKANTU_DEBUG_IN(); if (id == "potential") { Real epot = 0.; computePotentialEnergyByElements(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension), end = element_filter.lastType(spatial_dimension); for (; it != end ; ++it) { FEEngine & interface_engine = model->getFEEngine("EmbeddedInterfaceFEEngine"); epot += interface_engine.integrate(potential_energy(*it, _not_ghost), *it, _not_ghost, element_filter(*it, _not_ghost)); epot *= area; } return epot; } AKANTU_DEBUG_OUT(); return 0; } // Author is Nicolas Richart, see material.cc template void MaterialReinforcement::flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type, ElementKind element_kind) { AKANTU_DEBUG_IN(); typedef ElementTypeMapArray::type_iterator iterator; iterator tit = this->element_filter.firstType(1, ghost_type, element_kind); iterator end = this->element_filter.lastType(1, ghost_type, element_kind); for (; tit != end; ++tit) { ElementType type = *tit; try { __attribute__((unused)) const Array & src_vect = this->getArray(field_id,type,ghost_type); } catch(debug::Exception & e) { continue; } const Array & src_vect = this->getArray(field_id,type,ghost_type); const Array & filter = this->element_filter(type,ghost_type); // total number of elements for a given type UInt nb_element = this->model->getInterfaceMesh().getNbElement(type,ghost_type); // number of filtered elements UInt nb_element_src = filter.getSize(); // number of quadrature points per elem UInt nb_quad_per_elem = 0; // number of data per quadrature point UInt nb_data_per_quad = src_vect.getNbComponent(); if (!internal_flat.exists(type,ghost_type)) { internal_flat.alloc(nb_element*nb_quad_per_elem,nb_data_per_quad,type,ghost_type); } if (nb_element_src == 0) continue; nb_quad_per_elem = (src_vect.getSize()/nb_element_src); // number of data per element UInt nb_data = nb_quad_per_elem * src_vect.getNbComponent(); Array & dst_vect = internal_flat(type,ghost_type); dst_vect.resize(nb_element*nb_quad_per_elem); Array::const_scalar_iterator it = filter.begin(); Array::const_scalar_iterator end = filter.end(); Array::const_vector_iterator it_src = src_vect.begin_reinterpret(nb_data,nb_element_src); Array::vector_iterator it_dst = dst_vect.begin_reinterpret(nb_data,nb_element); for (; it != end ; ++it,++it_src) { it_dst[*it] = *it_src; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialReinforcement); /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh index ec22703b8..c7f9511b8 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh @@ -1,222 +1,228 @@ /** * @file material_reinforcement.hh * * @author Lucas Frérot * * @date creation: Thu Mar 12 2015 * @date last modification: Thu Mar 12 2015 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_REINFORCEMENT_HH__ #define __AKANTU_MATERIAL_REINFORCEMENT_HH__ #include "aka_common.hh" #include "material.hh" #include "embedded_interface_model.hh" #include "embedded_internal_field.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * @brief Material used to represent embedded reinforcements * * This class is used for computing the reinforcement stiffness matrix * along with the reinforcement residual. Room is made for constitutive law, * but actual use of contitutive laws is made in MaterialReinforcementTemplate. * * Be careful with the dimensions in this class : * - this->spatial_dimension is always 1 * - the template parameter dim is the dimension of the problem */ template class MaterialReinforcement : virtual public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Constructor MaterialReinforcement(SolidMechanicsModel & model, UInt spatial_dimension, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); /// Destructor virtual ~MaterialReinforcement(); protected: void initialize(SolidMechanicsModel & a_model); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Init the material virtual void initMaterial(); /// Init the background shape derivatives void initBackgroundShapeDerivatives(); /// Init the cosine matrices void initDirectingCosines(); /// Assemble stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// Update the residual virtual void updateResidual(GhostType ghost_type = _not_ghost); /// Assembled the residual virtual void assembleResidual(GhostType ghost_type); /// Compute all the stresses ! virtual void computeAllStresses(GhostType ghost_type); /// Compute the stiffness parameter for elements of a type virtual void computeTangentModuli(const ElementType & type, Array & tangent, GhostType ghost_type) = 0; + /// Compute energy virtual Real getEnergy(std::string id); + /// Flatten internal (modified from Material for 1D internals) void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type, ElementKind element_kind); /* ------------------------------------------------------------------------ */ /* Protected methods */ /* ------------------------------------------------------------------------ */ protected: - /// Allocate the background shape derivatives + /** + * @brief Allocate the background shape derivatives + * + * Background shape derivatives need to be stored per background element + * types but also per embedded element type, which is why they are stored + * in an ElementTypeMap *>. The outer ElementTypeMap + * refers to the embedded types, and the inner refers to the background types. + */ void allocBackgroundShapeDerivatives(); /// Compute the directing cosines matrix for one element type void computeDirectingCosines(const ElementType & type, GhostType ghost_type); /** * @brief Compute the directing cosines matrix on quadrature points. * * The structure of the directing cosines matrix is : * \f{eqnarray*}{ * C_{1,\cdot} & = & (l^2, m^2, n^2, lm, mn, ln) \\ * C_{i,j} & = & 0 * \f} * * with : * \f[ * (l, m, n) = \frac{1}{\|\frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}\|} \cdot \frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s} * \f] */ inline void computeDirectingCosinesOnQuad(const Matrix & nodes, Matrix & cosines); /// Assemble the stiffness matrix for an element type (typically _segment_2) void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /** * @brief Assemble the stiffness matrix for background & interface types * * Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001) * \f[ * \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T * \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}} * \f] */ void assembleStiffnessMatrix(const ElementType & interface_type, const ElementType & background_type, - GhostType interface_ghost, - GhostType background_ghost); + GhostType ghost_type); /// Compute the background shape derivatives for a type void computeBackgroundShapeDerivatives(const ElementType & type, GhostType ghost_type); /// Filter elements crossed by interface of a type void filterInterfaceBackgroundElements(Array & filter, const ElementType & type, const ElementType & interface_type, - GhostType ghost_type, - GhostType interface_ghost_type); + GhostType ghost_type); /// Assemble the residual of one type of element (typically _segment_2) void assembleResidual(const ElementType & type, GhostType ghost_type); /** * @brief Assemble the residual for a pair of elements * * Computes and assemble the residual. Residual in reinforcement is computed as : * \f[ * \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s} * \f] */ void assembleResidual(const ElementType & interface_type, const ElementType & background_type, - GhostType interface_ghost, - GhostType background_ghost); + GhostType ghost_type); // TODO figure out why voigt size is 4 in 2D inline void stressTensorToVoigtVector(const Matrix & tensor, Vector & vector); inline void strainTensorToVoigtVector(const Matrix & tensor, Vector & vector); /// Compute gradu on the interface quadrature points virtual void computeGradU(const ElementType & type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Embedded model EmbeddedInterfaceModel * model; /// Stress in the reinforcement InternalField stress_embedded; /// Gradu of concrete on reinforcement InternalField gradu_embedded; /// C matrix on quad InternalField directing_cosines; /// Prestress on quad InternalField pre_stress; /// Cross-sectional area Real area; /// Background mesh shape derivatives ElementTypeMap< ElementTypeMapArray * > shape_derivatives; }; #include "material_reinforcement_inline_impl.cc" __END_AKANTU__ #endif // __AKANTU_MATERIAL_REINFORCEMENT_HH__ diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc index a9ffebed4..8bde1600e 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_inline_impl.cc @@ -1,99 +1,99 @@ /** * @file material_reinforcement_inline_impl.cc * * @author Lucas Frérot * * @date creation: Thu Mar 12 2015 * @date last modification: Thu Mar 12 2015 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 2010-2012, 2014 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 . * */ /* -------------------------------------------------------------------------- */ -template -inline void MaterialReinforcement::computeDirectingCosinesOnQuad(const Matrix & nodes, +template +inline void MaterialReinforcement::computeDirectingCosinesOnQuad(const Matrix & nodes, Matrix & cosines) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(nodes.cols() == 2, "Higher order reinforcement elements not implemented"); const Vector & a = nodes(0), b = nodes(1); cosines.clear(); Real sq_length = 0.; - for (UInt i = 0 ; i < d ; i++) { + for (UInt i = 0 ; i < dim ; i++) { sq_length += Math::pow<2, Real>(b(i) - a(i)); } // Fill the first row of cosine matrix - for (UInt i = 0 ; i < d ; i++) { + for (UInt i = 0 ; i < dim ; i++) { cosines(0, i) = Math::pow<2, Real>(b(i) - a(i)) / sq_length; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialReinforcement::stressTensorToVoigtVector(const Matrix & tensor, Vector & vector) { AKANTU_DEBUG_IN(); for (UInt i = 0; i < dim; i++) { vector(i) = tensor(i, i); } if (dim == 2) { vector(2) = tensor(0, 1); } else if (dim == 3) { vector(3) = tensor(0, 1); vector(4) = tensor(0, 2); vector(5) = tensor(1, 2); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialReinforcement::strainTensorToVoigtVector(const Matrix & tensor, Vector & vector) { AKANTU_DEBUG_IN(); for (UInt i = 0; i < dim; i++) { vector(i) = tensor(i, i); } if (dim == 2) { vector(2) = 2 * tensor(0, 1); } else if (dim == 3) { vector(3) = 2 * tensor(0, 1); vector(4) = 2 * tensor(0, 2); vector(5) = 2 * tensor(1, 2); } AKANTU_DEBUG_OUT(); } diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_template_tmpl.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_template_tmpl.hh index 401f62f4f..e9ed3071f 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_template_tmpl.hh +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_template_tmpl.hh @@ -1,198 +1,197 @@ /** * @file material_reinforcement_template_tmpl.hh * * @author Lucas Frérot * * @date creation: Mon Mar 16 2015 * @date last modification: Mon Mar 16 2015 * * @brief Reinforcement material templated with constitutive law * * @section LICENSE * * Copyright (©) 2010-2012, 2014 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 . * */ // /!\ no namespace here ! /* -------------------------------------------------------------------------- */ template MaterialReinforcementTemplate::MaterialReinforcementTemplate(SolidMechanicsModel & a_model, const ID & id): Material(a_model, 1, dynamic_cast(a_model).getInterfaceMesh(), a_model.getFEEngine("EmbeddedInterfaceFEEngine"), id), MaterialReinforcement(a_model, 1, dynamic_cast(a_model).getInterfaceMesh(), a_model.getFEEngine("EmbeddedInterfaceFEEngine"), id), ConstLaw(a_model, 1, dynamic_cast(a_model).getInterfaceMesh(), a_model.getFEEngine("EmbeddedInterfaceFEEngine"), id) {} /* -------------------------------------------------------------------------- */ template MaterialReinforcementTemplate::~MaterialReinforcementTemplate() {} /* -------------------------------------------------------------------------- */ -/// TODO this function is super ugly. Find better way to initialize internals template void MaterialReinforcementTemplate::initMaterial() { MaterialReinforcement::initMaterial(); ConstLaw::initMaterial(); // Needed for plasticity law this->ConstLaw::nu = 0.5; this->ConstLaw::updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcementTemplate::computeGradU(const ElementType & el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialReinforcement::computeGradU(el_type, ghost_type); const UInt voigt_size = Material::getTangentStiffnessVoigtSize(dim); Array::matrix_iterator full_gradu_it = this->MaterialReinforcement::gradu_embedded(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator full_gradu_end = this->MaterialReinforcement::gradu_embedded(el_type, ghost_type).end(dim, dim); Array::scalar_iterator gradu_it = this->ConstLaw::gradu(el_type, ghost_type).begin(); Array::matrix_iterator cosines_it = this->directing_cosines(el_type, ghost_type).begin(voigt_size, voigt_size); for (; full_gradu_it != full_gradu_end ; ++full_gradu_it, ++gradu_it, ++cosines_it) { Matrix & full_gradu = *full_gradu_it; Real & gradu = *gradu_it; Matrix & C = *cosines_it; computeInterfaceGradUOnQuad(full_gradu, gradu, C); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcementTemplate::computeInterfaceGradUOnQuad(const Matrix & full_gradu, Real & gradu, const Matrix & C) { const UInt voigt_size = Material::getTangentStiffnessVoigtSize(dim); Matrix epsilon(dim, dim); Vector e_voigt(voigt_size); Vector e_interface_voigt(voigt_size); epsilon = 0.5 * (full_gradu + full_gradu.transpose()); MaterialReinforcement::strainTensorToVoigtVector(epsilon, e_voigt); e_interface_voigt.mul(C, e_voigt); gradu = e_interface_voigt(0); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcementTemplate::computeTangentModuli(const ElementType & el_type, Array & tangent, GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(tangent.getNbComponent() == 1, "Reinforcements only work in 1D"); ConstLaw::computeTangentModuli(el_type, tangent, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcementTemplate::computeStress(ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); ConstLaw::computeStress(type, ghost_type); Array::matrix_iterator full_sigma_it = this->MaterialReinforcement::stress_embedded(type, ghost_type).begin(dim, dim); Array::matrix_iterator full_sigma_end = this->MaterialReinforcement::stress_embedded(type, ghost_type).end(dim, dim); Array::scalar_iterator sigma_it = this->ConstLaw::stress(type, ghost_type).begin(); Array::scalar_iterator pre_stress_it = this->MaterialReinforcement::pre_stress(type, ghost_type).begin(); for (; full_sigma_it != full_sigma_end ; ++full_sigma_it, ++sigma_it, ++pre_stress_it) { Matrix & sigma = *full_sigma_it; sigma(0, 0) = *sigma_it + *pre_stress_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real MaterialReinforcementTemplate::getEnergy(std::string id) { return MaterialReinforcement::getEnergy(id); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcementTemplate::computePotentialEnergy(ElementType type, GhostType ghost_type) { const UInt nb_elements = this->element_filter(type, ghost_type).getSize(); const UInt nb_quad = this->model->getFEEngine("EmbeddedInterfaceFEEngine").getNbQuadraturePoints(type); this->ConstLaw::potential_energy.alloc(nb_quad * nb_elements, 1, type, ghost_type, 0.); ConstLaw::computePotentialEnergy(type, ghost_type); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcementTemplate::flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type, ElementKind element_kind) { if (field_id == "stress_embedded" || field_id == "inelastic_strain") MaterialReinforcement::flattenInternal(field_id, internal_flat, ghost_type, element_kind); } template void MaterialReinforcementTemplate::savePreviousState() { ConstLaw::savePreviousState(); }