diff --git a/python/swig/mesh.i b/python/swig/mesh.i index eb036f68d..f294f8730 100644 --- a/python/swig/mesh.i +++ b/python/swig/mesh.i @@ -1,231 +1,230 @@ /** * @file mesh.i * * @author Guillaume Anciaux * @author Fabian Barras * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Dec 12 2014 * @date last modification: Wed Jan 13 2016 * * @brief mesh wrapper * * @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 . * */ %{ #include "mesh.hh" #include "node_group.hh" #include "solid_mechanics_model.hh" #include "python_functor.hh" #include "mesh_utils.hh" #include "aka_bbox.hh" #include "mesh_accessor.hh" #include "communicator.hh" using akantu::IntegrationPoint; using akantu::Vector; using akantu::ElementTypeMapArray; using akantu::MatrixProxy; using akantu::Matrix; using akantu::UInt; using akantu::Real; using akantu::Array; using akantu::BBox; using akantu::Communicator; using akantu::SolidMechanicsModel; %} namespace akantu { %ignore NewNodesEvent; %ignore RemovedNodesEvent; %ignore NewElementsEvent; %ignore RemovedElementsEvent; %ignore MeshEventHandler; %ignore MeshEvent< UInt >; %ignore MeshEvent< Element >; %ignore Mesh::extractNodalCoordinatesFromPBCElement; %ignore Mesh::getGroupDumer; %ignore Mesh::getFacetLocalConnectivity; %ignore Mesh::getAllFacetTypes; %ignore Mesh::getCommunicator; %ignore Mesh::getConnectivities; %ignode Mesh::getBBox; %ignore GroupManager::getElementGroups; %ignore Dumpable::addDumpFieldExternalReal; } print_self(Mesh) // Swig considers enums to be ints, and it creates a conflict with two versions of getNbElement() %rename(getNbElementByDimension) akantu::Mesh::getNbElement(const UInt spatial_dimension = _all_dimensions, const GhostType& ghost_type = _not_ghost, const ElementKind& kind = _ek_not_defined) const; %extend akantu::Mesh { - - PyObject * getElementGroups(){ - return akantu::PythonFunctor::convertToPython($self->getElementGroups()); - } + // PyObject * getElementGroups(){ + // return akantu::PythonFunctor::convertToPython($self->getElementGroups()); + // } PyObject * getAllConnectivities(){ return akantu::PythonFunctor::convertToPython($self->getConnectivities()); } void resizeMesh(UInt nb_nodes, UInt nb_element, const ElementType & type) { Array & nodes = const_cast &>($self->getNodes()); nodes.resize(nb_nodes); $self->addConnectivityType(type); Array & connectivity = const_cast &>($self->getConnectivity(type)); connectivity.resize(nb_element); } Array & getNodalDataReal(const ID & name, UInt nb_components = 1) { auto && data = $self->getNodalData(name, nb_components); data.resize($self->getNbNodes()); return data; } bool hasDataReal(const ID & name, const ElementType & type) { return $self->hasData(name, type); } Array & getElementalDataReal(const ID & name, const ElementType & type, UInt nb_components = 1) { auto && data = $self->getElementalDataArrayAlloc(name, type, akantu::_not_ghost, nb_components); data.resize($self->getNbElement(type, akantu::_not_ghost)); return data; } Array & getElementalDataUInt(const ID & name, const ElementType & type, UInt nb_components = 1) { auto && data = $self->getElementalDataArrayAlloc(name, type, akantu::_not_ghost, nb_components); data.resize($self->getNbElement(type, akantu::_not_ghost)); return data; } Array & computeBarycenters(const ElementType & type) { auto dim = $self->getSpatialDimension(); auto && data = $self->getElementalDataArrayAlloc("barycenters", type, akantu::_not_ghost, dim); auto nb_el = data.size(); auto total_nb_el = $self->getNbElement(type, akantu::_not_ghost); data.resize(total_nb_el); auto bary_it = make_view(data, dim).begin() + nb_el; for (auto el = nb_el; el < total_nb_el; ++el) { $self->getBarycenter(akantu::Element{type, el, akantu::_not_ghost}, *bary_it); ++bary_it; } return data; } void ready() { akantu::MeshAccessor ma(* $self); ma.makeReady(); } } %extend akantu::GroupManager { void createGroupsFromStringMeshData(const std::string & dataset_name) { if (dataset_name == "physical_names"){ AKANTU_EXCEPTION("Deprecated behavior: no need to call 'createGroupsFromStringMeshData' for physical names"); } $self->createGroupsFromMeshData(dataset_name); } void createGroupsFromUIntMeshData(const std::string & dataset_name) { $self->createGroupsFromMeshData(dataset_name); } } %extend akantu::NodeGroup { akantu::Array & getGroupedNodes(akantu::Array & surface_array, Mesh & mesh) { auto && group_node = $self->getNodes(); auto && full_array = mesh.getNodes(); surface_array.resize(group_node.size()); for (UInt i = 0; i < group_node.size(); ++i) { for (UInt cmp = 0; cmp < full_array.getNbComponent(); ++cmp) { surface_array(i, cmp) = full_array(group_node(i), cmp); } } akantu::Array & res(surface_array); return res; } akantu::Array & getGroupedArray(akantu::Array & surface_array, akantu::SolidMechanicsModel & model, int type) { akantu::Array * full_array; switch (type) { case 0 : full_array = new akantu::Array(model.getDisplacement()); break; case 1 : full_array = new akantu::Array(model.getVelocity()); break; case 2 : full_array = new akantu::Array(model.getExternalForce()); break; } akantu::Array group_node = $self->getNodes(); surface_array.resize(group_node.size()); for (UInt i = 0; i < group_node.size(); ++i) { for (UInt cmp = 0; cmp < full_array->getNbComponent(); ++cmp) { surface_array(i,cmp) = (*full_array)(group_node(i),cmp); } } akantu::Array & res(surface_array); return res; } } %include "group_manager.hh" %include "node_group.hh" %include "dumper_iohelper.hh" %include "dumpable_iohelper.hh" %include "element_group.hh" %include "mesh.hh" %include "mesh_utils.hh" %include "aka_bbox.hh" namespace akantu{ %extend Dumpable { void addDumpFieldExternalReal(const std::string & field_id, const Array & field){ $self->addDumpFieldExternal(field_id,field); } } } diff --git a/src/io/mesh_io/mesh_io_msh.cc b/src/io/mesh_io/mesh_io_msh.cc index b1781bee4..57102de4f 100644 --- a/src/io/mesh_io/mesh_io_msh.cc +++ b/src/io/mesh_io/mesh_io_msh.cc @@ -1,1069 +1,1089 @@ /** * @file mesh_io_msh.cc * * @author Dana Christen * @author Mauro Corrado * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief Read/Write for MSH files generated by gmsh * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* ----------------------------------------------------------------------------- Version (Legacy) 1.0 $NOD number-of-nodes node-number x-coord y-coord z-coord ... $ENDNOD $ELM number-of-elements elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list ... $ENDELM ----------------------------------------------------------------------------- Version 2.1 $MeshFormat version-number file-type data-size $EndMeshFormat $Nodes number-of-nodes node-number x-coord y-coord z-coord ... $EndNodes $Elements number-of-elements elm-number elm-type number-of-tags < tag > ... node-number-list ... $EndElements $PhysicalNames number-of-names physical-dimension physical-number "physical-name" ... $EndPhysicalNames $NodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... node-number value ... ... $EndNodeData $ElementData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number value ... ... $EndElementData $ElementNodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number number-of-nodes-per-element value ... ... $ElementEndNodeData ----------------------------------------------------------------------------- elem-type 1: 2-node line. 2: 3-node triangle. 3: 4-node quadrangle. 4: 4-node tetrahedron. 5: 8-node hexahedron. 6: 6-node prism. 7: 5-node pyramid. 8: 3-node second order line 9: 6-node second order triangle 10: 9-node second order quadrangle 11: 10-node second order tetrahedron 12: 27-node second order hexahedron 13: 18-node second order prism 14: 14-node second order pyramid 15: 1-node point. 16: 8-node second order quadrangle 17: 20-node second order hexahedron 18: 15-node second order prism 19: 13-node second order pyramid 20: 9-node third order incomplete triangle 21: 10-node third order triangle 22: 12-node fourth order incomplete triangle 23: 15-node fourth order triangle 24: 15-node fifth order incomplete triangle 25: 21-node fifth order complete triangle 26: 4-node third order edge 27: 5-node fourth order edge 28: 6-node fifth order edge 29: 20-node third order tetrahedron 30: 35-node fourth order tetrahedron 31: 56-node fifth order tetrahedron -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "mesh_io.hh" #include "mesh_utils.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Methods Implentations */ /* -------------------------------------------------------------------------- */ MeshIOMSH::MeshIOMSH() { canReadSurface = true; canReadExtendedData = true; _msh_nodes_per_elem[_msh_not_defined] = 0; _msh_nodes_per_elem[_msh_segment_2] = 2; _msh_nodes_per_elem[_msh_triangle_3] = 3; _msh_nodes_per_elem[_msh_quadrangle_4] = 4; _msh_nodes_per_elem[_msh_tetrahedron_4] = 4; _msh_nodes_per_elem[_msh_hexahedron_8] = 8; _msh_nodes_per_elem[_msh_prism_1] = 6; _msh_nodes_per_elem[_msh_pyramid_1] = 1; _msh_nodes_per_elem[_msh_segment_3] = 3; _msh_nodes_per_elem[_msh_triangle_6] = 6; _msh_nodes_per_elem[_msh_quadrangle_9] = 9; _msh_nodes_per_elem[_msh_tetrahedron_10] = 10; _msh_nodes_per_elem[_msh_hexahedron_27] = 27; _msh_nodes_per_elem[_msh_hexahedron_20] = 20; _msh_nodes_per_elem[_msh_prism_18] = 18; _msh_nodes_per_elem[_msh_prism_15] = 15; _msh_nodes_per_elem[_msh_pyramid_14] = 14; _msh_nodes_per_elem[_msh_point] = 1; _msh_nodes_per_elem[_msh_quadrangle_8] = 8; _msh_to_akantu_element_types[_msh_not_defined] = _not_defined; _msh_to_akantu_element_types[_msh_segment_2] = _segment_2; _msh_to_akantu_element_types[_msh_triangle_3] = _triangle_3; _msh_to_akantu_element_types[_msh_quadrangle_4] = _quadrangle_4; _msh_to_akantu_element_types[_msh_tetrahedron_4] = _tetrahedron_4; _msh_to_akantu_element_types[_msh_hexahedron_8] = _hexahedron_8; _msh_to_akantu_element_types[_msh_prism_1] = _pentahedron_6; _msh_to_akantu_element_types[_msh_pyramid_1] = _not_defined; _msh_to_akantu_element_types[_msh_segment_3] = _segment_3; _msh_to_akantu_element_types[_msh_triangle_6] = _triangle_6; _msh_to_akantu_element_types[_msh_quadrangle_9] = _not_defined; _msh_to_akantu_element_types[_msh_tetrahedron_10] = _tetrahedron_10; _msh_to_akantu_element_types[_msh_hexahedron_27] = _not_defined; _msh_to_akantu_element_types[_msh_hexahedron_20] = _hexahedron_20; _msh_to_akantu_element_types[_msh_prism_18] = _not_defined; _msh_to_akantu_element_types[_msh_prism_15] = _pentahedron_15; _msh_to_akantu_element_types[_msh_pyramid_14] = _not_defined; _msh_to_akantu_element_types[_msh_point] = _point_1; _msh_to_akantu_element_types[_msh_quadrangle_8] = _quadrangle_8; _akantu_to_msh_element_types[_not_defined] = _msh_not_defined; _akantu_to_msh_element_types[_segment_2] = _msh_segment_2; _akantu_to_msh_element_types[_segment_3] = _msh_segment_3; _akantu_to_msh_element_types[_triangle_3] = _msh_triangle_3; _akantu_to_msh_element_types[_triangle_6] = _msh_triangle_6; _akantu_to_msh_element_types[_tetrahedron_4] = _msh_tetrahedron_4; _akantu_to_msh_element_types[_tetrahedron_10] = _msh_tetrahedron_10; _akantu_to_msh_element_types[_quadrangle_4] = _msh_quadrangle_4; _akantu_to_msh_element_types[_quadrangle_8] = _msh_quadrangle_8; _akantu_to_msh_element_types[_hexahedron_8] = _msh_hexahedron_8; _akantu_to_msh_element_types[_hexahedron_20] = _msh_hexahedron_20; _akantu_to_msh_element_types[_pentahedron_6] = _msh_prism_1; _akantu_to_msh_element_types[_pentahedron_15] = _msh_prism_15; _akantu_to_msh_element_types[_point_1] = _msh_point; #if defined(AKANTU_STRUCTURAL_MECHANICS) _akantu_to_msh_element_types[_bernoulli_beam_2] = _msh_segment_2; _akantu_to_msh_element_types[_bernoulli_beam_3] = _msh_segment_2; _akantu_to_msh_element_types[_discrete_kirchhoff_triangle_18] = _msh_triangle_3; #endif std::map::iterator it; for (it = _akantu_to_msh_element_types.begin(); it != _akantu_to_msh_element_types.end(); ++it) { UInt nb_nodes = _msh_nodes_per_elem[it->second]; std::vector tmp(nb_nodes); for (UInt i = 0; i < nb_nodes; ++i) { tmp[i] = i; } switch (it->first) { case _tetrahedron_10: tmp[8] = 9; tmp[9] = 8; break; case _pentahedron_6: tmp[0] = 2; tmp[1] = 0; tmp[2] = 1; tmp[3] = 5; tmp[4] = 3; tmp[5] = 4; break; case _pentahedron_15: tmp[0] = 2; tmp[1] = 0; tmp[2] = 1; tmp[3] = 5; tmp[4] = 3; tmp[5] = 4; tmp[6] = 8; tmp[8] = 11; tmp[9] = 6; tmp[10] = 9; tmp[11] = 10; tmp[12] = 14; tmp[14] = 12; break; case _hexahedron_20: tmp[9] = 11; tmp[10] = 12; tmp[11] = 9; tmp[12] = 13; tmp[13] = 10; tmp[17] = 19; tmp[18] = 17; tmp[19] = 18; break; default: // nothing to change break; } _read_order[it->first] = tmp; } } /* -------------------------------------------------------------------------- */ MeshIOMSH::~MeshIOMSH() = default; /* -------------------------------------------------------------------------- */ namespace { struct File { std::string filename; std::ifstream infile; std::string line; size_t current_line{0}; size_t first_node_number{std::numeric_limits::max()}, last_node_number{0}; bool has_physical_names{false}; std::unordered_map node_tags; std::unordered_map element_tags; double version{0}; int size_of_size_t{0}; Mesh & mesh; MeshAccessor mesh_accessor; - std::multimap entity_tag_to_physical_tags; + std::multimap, int> entity_tag_to_physical_tags; File(const std::string & filename, Mesh & mesh) : filename(filename), mesh(mesh), mesh_accessor(mesh) { infile.open(filename.c_str()); if (not infile.good()) { AKANTU_EXCEPTION("Cannot open file " << filename); } } ~File() { infile.close(); } auto good() { return infile.good(); } std::stringstream get_line() { std::string tmp_str; if (infile.eof()) { AKANTU_EXCEPTION("Reached the end of the file " << filename); } std::getline(infile, tmp_str); line = trim(tmp_str); ++current_line; return std::stringstream(line); } template void read_line(Ts &&... ts) { auto && sstr = get_line(); (void)std::initializer_list{ (sstr >> std::forward(ts), 0)...}; } }; } // namespace /* -------------------------------------------------------------------------- */ template void MeshIOMSH::populateReaders2(File & file, Readers & readers) { readers["$NOD"] = readers["$Nodes"] = [&](const std::string &) { UInt nb_nodes; file.read_line(nb_nodes); Array & nodes = file.mesh_accessor.getNodes(); nodes.resize(nb_nodes); file.mesh_accessor.setNbGlobalNodes(nb_nodes); size_t index; Vector coord(3); /// for each node, read the coordinates for (auto && data : enumerate(make_view(nodes, nodes.getNbComponent()))) { file.read_line(index, coord(0), coord(1), coord(2)); if (index > std::numeric_limits::max()) { AKANTU_EXCEPTION( "There are more nodes in this files than the index type of akantu " "can handle, consider recompiling with a bigger index type"); } file.first_node_number = std::min(file.first_node_number, index); file.last_node_number = std::max(file.last_node_number, index); for (auto && coord_data : zip(std::get<1>(data), coord)) { std::get<0>(coord_data) = std::get<1>(coord_data); } file.node_tags[index] = std::get<0>(data); } }; readers["$ELM"] = readers["$Elements"] = [&](const std::string &) { UInt nb_elements; file.read_line(nb_elements); Int index; UInt msh_type; ElementType akantu_type; for (UInt i = 0; i < nb_elements; ++i) { auto && sstr_elem = file.get_line(); sstr_elem >> index; sstr_elem >> msh_type; /// get the connectivity vector depending on the element type akantu_type = this->_msh_to_akantu_element_types[MSHElementType(msh_type)]; if (akantu_type == _not_defined) { AKANTU_DEBUG_WARNING("Unsuported element kind " << msh_type << " at line " << file.current_line); continue; } Element elem{akantu_type, 0, _not_ghost}; auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type); auto node_per_element = connectivity.getNbComponent(); auto & read_order = this->_read_order[akantu_type]; /// read tags informations if (file.version < 2) { Int tag0, tag1, nb_nodes; // reg-phys, reg-elem, number-of-nodes sstr_elem >> tag0 >> tag1 >> nb_nodes; auto & data0 = file.mesh_accessor.template getData("tag_0", akantu_type); data0.push_back(tag0); auto & data1 = file.mesh_accessor.template getData("tag_1", akantu_type); data1.push_back(tag1); } else if (file.version < 4) { UInt nb_tags; sstr_elem >> nb_tags; for (UInt j = 0; j < nb_tags; ++j) { Int tag; sstr_elem >> tag; auto & data = file.mesh_accessor.template getData( "tag_" + std::to_string(j), akantu_type); data.push_back(tag); } } Vector local_connect(node_per_element); for (UInt j = 0; j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; AKANTU_DEBUG_ASSERT(node_index <= file.last_node_number, "Node number not in range : line " << file.current_line); local_connect(read_order[j]) = file.node_tags[node_index]; } connectivity.push_back(local_connect); elem.element = connectivity.size() - 1; file.element_tags[index] = elem; } }; readers["$Periodic"] = [&](const std::string &) { UInt nb_periodic_entities; file.read_line(nb_periodic_entities); file.mesh_accessor.getNodesFlags().resize(file.mesh.getNbNodes(), NodeFlag::_normal); for (UInt p = 0; p < nb_periodic_entities; ++p) { // dimension slave-tag master-tag UInt dimension; file.read_line(dimension); // transformation file.get_line(); // nb nodes UInt nb_nodes; file.read_line(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { // slave master auto && sstr = file.get_line(); // The info in the mesh seem inconsistent so they are ignored for now. continue; if (dimension == file.mesh.getSpatialDimension() - 1) { UInt slave, master; sstr >> slave; sstr >> master; file.mesh_accessor.addPeriodicSlave(file.node_tags[slave], file.node_tags[master]); } } } // mesh_accessor.markMeshPeriodic(); }; } /* -------------------------------------------------------------------------- */ template void MeshIOMSH::populateReaders4(File & file, Readers & readers) { static std::map entity_type{ {0, "points"}, {1, "curve"}, {2, "surface"}, {3, "volume"}, }; readers["$Entities"] = [&](const std::string &) { size_t num_entity[4]; file.read_line(num_entity[0], num_entity[1], num_entity[2], num_entity[3]); for (auto entity_dim : arange(4)) { for (auto _ [[gnu::unused]] : arange(num_entity[entity_dim])) { auto && sstr = file.get_line(); int tag; double min_x, min_y, min_z, max_x, max_y, max_z; size_t num_physical_tags; sstr >> tag >> min_x >> min_y >> min_z; if (entity_dim > 0 or file.version < 4.1) { sstr >> max_x >> max_y >> max_z; } sstr >> num_physical_tags; for (auto _ [[gnu::unused]] : arange(num_physical_tags)) { int phys_tag; sstr >> phys_tag; std::string physical_name; if (this->physical_names.find(phys_tag) == this->physical_names.end()) { physical_name = "msh_block_" + std::to_string(phys_tag); } else { physical_name = this->physical_names[phys_tag]; } if (not file.mesh.elementGroupExists(physical_name)) { file.mesh.createElementGroup(physical_name, entity_dim); } else { file.mesh.getElementGroup(physical_name).addDimension(entity_dim); } file.entity_tag_to_physical_tags.insert( - std::make_pair(tag, phys_tag)); + std::make_pair(std::make_pair(tag, entity_dim), phys_tag)); } } } }; readers["$Nodes"] = [&](const std::string &) { size_t num_blocks, num_nodes; if (file.version >= 4.1) { file.read_line(num_blocks, num_nodes, file.first_node_number, file.last_node_number); } else { file.read_line(num_blocks, num_nodes); } auto & nodes = file.mesh_accessor.getNodes(); nodes.reserve(num_nodes); file.mesh_accessor.setNbGlobalNodes(num_nodes); if (num_nodes > std::numeric_limits::max()) { AKANTU_EXCEPTION( "There are more nodes in this files than the index type of akantu " "can handle, consider recompiling with a bigger index type"); } + size_t node_id{0}; + for (auto block [[gnu::unused]] : arange(num_blocks)) { int entity_dim, entity_tag, parametric; size_t num_nodes_in_block; - size_t node_id{0}; - - // auto & grp = file.mesh.createNodeGroup("msh_" + entity_type[entity_dim] - // + - // "_block_" + - // std::to_string(block)); + Vector pos(3); + Vector real_pos(nodes.getNbComponent()); if (file.version >= 4.1) { file.read_line(entity_dim, entity_tag, parametric, num_nodes_in_block); if (parametric) { AKANTU_EXCEPTION( "Akantu does not support parametric nodes in msh files"); } for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) { size_t tag; file.read_line(tag); file.node_tags[tag] = node_id; - // grp.add(node_id); ++node_id; } for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) { - Vector pos(3); file.read_line(pos(_x), pos(_y), pos(_z)); - nodes.push_back(pos); + for (auto data : zip(real_pos, pos)) { + std::get<0>(data) = std::get<1>(data); + } + nodes.push_back(real_pos); } } else { file.read_line(entity_tag, entity_dim, parametric, num_nodes_in_block); for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) { size_t tag; - Vector pos(3); file.read_line(tag, pos(_x), pos(_y), pos(_z)); if (file.version < 4.1) { file.first_node_number = std::min(file.first_node_number, tag); file.last_node_number = std::max(file.last_node_number, tag); } - nodes.push_back(pos); + for (auto data : zip(real_pos, pos)) { + std::get<0>(data) = std::get<1>(data); + } + nodes.push_back(real_pos); file.node_tags[tag] = node_id; - // grp.add(node_id); ++node_id; } } } }; readers["$Elements"] = [&](const std::string &) { size_t num_blocks, num_elements; file.read_line(num_blocks, num_elements); for (auto block [[gnu::unused]] : arange(num_blocks)) { int entity_dim, entity_tag, element_type; size_t num_elements_in_block; if (file.version >= 4.1) { file.read_line(entity_dim, entity_tag, element_type, num_elements_in_block); } else { file.read_line(entity_tag, entity_dim, element_type, num_elements_in_block); } /// get the connectivity vector depending on the element type auto && akantu_type = this->_msh_to_akantu_element_types[(MSHElementType)element_type]; if (akantu_type == _not_defined) { AKANTU_DEBUG_WARNING("Unsuported element kind " << element_type << " at line " << file.current_line); continue; } Element elem{akantu_type, 0, _not_ghost}; auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type); Vector local_connect(connectivity.getNbComponent()); auto && read_order = this->_read_order[akantu_type]; + auto & data0 = + file.mesh_accessor.template getData("tag_0", akantu_type); + data0.resize(data0.size() + num_elements_in_block, 0); + + auto & physical_data = file.mesh_accessor.template getData( + "physical_names", akantu_type); + physical_data.resize(physical_data.size() + num_elements_in_block, ""); + for (auto _ [[gnu::unused]] : arange(num_elements_in_block)) { auto && sstr_elem = file.get_line(); size_t elem_tag; sstr_elem >> elem_tag; for (auto && c : arange(connectivity.getNbComponent())) { size_t node_tag; sstr_elem >> node_tag; AKANTU_DEBUG_ASSERT(node_tag <= file.last_node_number, "Node number not in range : line " << file.current_line); node_tag = file.node_tags[node_tag]; local_connect(read_order[c]) = node_tag; } connectivity.push_back(local_connect); elem.element = connectivity.size() - 1; file.element_tags[elem_tag] = elem; - auto range = file.entity_tag_to_physical_tags.equal_range(entity_tag); + auto range = file.entity_tag_to_physical_tags.equal_range( + std::make_pair(entity_tag, entity_dim)); + bool first = true; for (auto it = range.first; it != range.second; ++it) { - file.mesh.getElementGroup(this->physical_names[it->second]).add(elem); + const auto & str = this->physical_names[it->second]; + if (first) { + data0(elem.element) = + it->second; // for compatibility with version 2 + physical_data(elem.element) = str; + first = false; + } + file.mesh.getElementGroup(str).add(elem, true, false); } } } + + for (auto && element_group : file.mesh.iterateElementGroups()) { + element_group.getNodeGroup().optimize(); + } }; } /* -------------------------------------------------------------------------- */ void MeshIOMSH::read(const std::string & filename, Mesh & mesh) { File file(filename, mesh); std::map> readers; readers["$MeshFormat"] = [&](const std::string &) { auto && sstr = file.get_line(); int format; sstr >> file.version >> format; if (format != 0) AKANTU_ERROR("This reader can only read ASCII files."); if (file.version > 2) { sstr >> file.size_of_size_t; if (file.size_of_size_t > int(sizeof(UInt))) { AKANTU_DEBUG_INFO("The size of the indexes in akantu might be to small " "to read this file (akantu " << sizeof(UInt) << " vs. msh file " << file.size_of_size_t << ")"); } } if (file.version < 4) { this->populateReaders2(file, readers); } else { this->populateReaders4(file, readers); } }; auto && read_data = [&](auto && entity_tags, auto && get_data, auto && read_data) { auto read_data_tags = [&](auto x) { UInt number_of_tags{0}; file.read_line(number_of_tags); std::vector tags(number_of_tags); for (auto && tag : tags) { file.read_line(tag); } return tags; }; auto && string_tags = read_data_tags(std::string{}); auto && real_tags [[gnu::unused]] = read_data_tags(double{}); auto && int_tags = read_data_tags(int{}); for (auto & s : string_tags) { s = trim(s, '"'); } auto id = string_tags[0]; auto size = int_tags[2]; auto nb_component = int_tags[1]; auto & data = get_data(id, size, nb_component); for (auto n [[gnu::unused]] : arange(size)) { auto && sstr = file.get_line(); size_t tag; sstr >> tag; const auto & entity = entity_tags[tag]; read_data(entity, sstr, data, nb_component); } }; readers["$NodeData"] = [&](const std::string &) { /* $NodeData numStringTags(ASCII int) stringTag(string) ... numRealTags(ASCII int) realTag(ASCII double) ... numIntegerTags(ASCII int) integerTag(ASCII int) ... nodeTag(size_t) value(double) ... ... $EndNodeData */ read_data(file.node_tags, [&](auto && id, auto && size [[gnu::unused]], auto && nb_component [[gnu::unused]]) -> Array & { auto & data = file.mesh.template registerNodalData( id, nb_component); data.resize(size); return data; }, [&](auto && node, auto && sstr, auto && data, auto && nb_component [[gnu::unused]]) { for (auto c : arange(nb_component)) { sstr >> data(node, c); } }); }; readers["$ElementData"] = [&](const std::string &) { /* $ElementData numStringTags(ASCII int) stringTag(string) ... numRealTags(ASCII int) realTag(ASCII double) ... numIntegerTags(ASCII int) integerTag(ASCII int) ... elementTag(size_t) value(double) ... ... $EndElementData */ read_data( file.element_tags, [&](auto && id, auto && size [[gnu::unused]], auto && nb_component [[gnu::unused]]) -> ElementTypeMapArray & { file.mesh.template registerElementalData(id); return file.mesh.template getElementalData(id); }, [&](auto && element, auto && sstr, auto && data, auto && nb_component) { if (not data.exists(element.type)) { data.alloc(mesh.getNbElement(element.type), nb_component, element.type, element.ghost_type); } auto & data_array = data(element.type); for (auto c : arange(nb_component)) { sstr >> data_array(element.element, c); } }); }; readers["$ElementNodeData"] = [&](const std::string &) { /* $ElementNodeData numStringTags(ASCII int) stringTag(string) ... numRealTags(ASCII int) realTag(ASCII double) ... numIntegerTags(ASCII int) integerTag(ASCII int) ... elementTag(size_t) value(double) ... ... $EndElementNodeData */ read_data( file.element_tags, [&](auto && id, auto && size [[gnu::unused]], auto && nb_component [[gnu::unused]]) -> ElementTypeMapArray & { file.mesh.template registerElementalData(id); auto & data = file.mesh.template getElementalData(id); data.isNodal(true); return data; }, [&](auto && element, auto && sstr, auto && data, auto && nb_component) { int nb_nodes_per_element; sstr >> nb_nodes_per_element; if (not data.exists(element.type)) { data.alloc(mesh.getNbElement(element.type), nb_component * nb_nodes_per_element, element.type, element.ghost_type); } auto & data_array = data(element.type); for (auto c : arange(nb_component)) { sstr >> data_array(element.element, c); } }); }; readers["$PhysicalNames"] = [&](const std::string &) { file.has_physical_names = true; int num_of_phys_names; file.read_line(num_of_phys_names); /// the format line for (auto k [[gnu::unused]] : arange(num_of_phys_names)) { int phys_name_id; int phys_dim; std::string phys_name; file.read_line(phys_dim, phys_name_id, std::quoted(phys_name)); this->physical_names[phys_name_id] = phys_name; } }; readers["Unsupported"] = [&](const std::string & _block) { std::string block = _block.substr(1); AKANTU_DEBUG_WARNING("Unsupported block_kind " << block << " at line " << file.current_line); auto && end_block = "$End" + block; while (file.line != end_block) { file.get_line(); } }; while (file.good()) { std::string block; file.read_line(block); auto && it = readers.find(block); if (it != readers.end()) { it->second(block); std::string end_block; file.read_line(end_block); block = block.substr(1); if (end_block != "$End" + block) { AKANTU_EXCEPTION("The reader failed to properly read the block " << block << ". Expected a $End" << block << " at line " << file.current_line); } } else if (block.size() != 0) { readers["Unsupported"](block); } } // mesh.updateTypesOffsets(_not_ghost); - if (file.version < 4) { this->constructPhysicalNames("tag_0", mesh); if (file.has_physical_names) mesh.createGroupsFromMeshData("physical_names"); } MeshUtils::fillElementToSubElementsData(mesh); } /* -------------------------------------------------------------------------- */ void MeshIOMSH::write(const std::string & filename, const Mesh & mesh) { std::ofstream outfile; const Array & nodes = mesh.getNodes(); outfile.open(filename.c_str()); outfile << "$MeshFormat" << "\n"; outfile << "2.2 0 8" << "\n"; outfile << "$EndMeshFormat" << "\n"; outfile << std::setprecision(std::numeric_limits::digits10); outfile << "$Nodes" << "\n"; outfile << nodes.size() << "\n"; outfile << std::uppercase; for (UInt i = 0; i < nodes.size(); ++i) { Int offset = i * nodes.getNbComponent(); outfile << i + 1; for (UInt j = 0; j < nodes.getNbComponent(); ++j) { outfile << " " << nodes.storage()[offset + j]; } for (UInt p = nodes.getNbComponent(); p < 3; ++p) outfile << " " << 0.; outfile << "\n"; ; } outfile << std::nouppercase; outfile << "$EndNodes" << "\n"; outfile << "$Elements" << "\n"; Int nb_elements = 0; for (auto && type : mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { const Array & connectivity = mesh.getConnectivity(type, _not_ghost); nb_elements += connectivity.size(); } outfile << nb_elements << "\n"; std::map element_to_msh_element; UInt element_idx = 1; Element element; for (auto && type : mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { const auto & connectivity = mesh.getConnectivity(type, _not_ghost); element.type = type; UInt * tag[2] = {nullptr, nullptr}; if (mesh.hasData("tag_0", type, _not_ghost)) { const auto & data_tag_0 = mesh.getData("tag_0", type, _not_ghost); tag[0] = data_tag_0.storage(); } if (mesh.hasData("tag_1", type, _not_ghost)) { const auto & data_tag_1 = mesh.getData("tag_1", type, _not_ghost); tag[1] = data_tag_1.storage(); } for (auto && data : enumerate(make_view(connectivity, connectivity.getNbComponent()))) { element.element = std::get<0>(data); const auto & conn = std::get<1>(data); element_to_msh_element[element] = element_idx; outfile << element_idx << " " << _akantu_to_msh_element_types[type] << " 2"; /// \todo write the real data in the file for (UInt t = 0; t < 2; ++t) if (tag[t]) outfile << " " << tag[t][element.element]; else outfile << " 0"; for (auto && c : conn) { outfile << " " << c + 1; } outfile << "\n"; element_idx++; } } outfile << "$EndElements" << "\n"; if (mesh.hasData(MeshDataType::_nodal)) { auto && tags = mesh.getTagNames(); for (auto && tag : tags) { auto type = mesh.getTypeCode(tag, MeshDataType::_nodal); if (type != MeshDataTypeCode::_real) { AKANTU_DEBUG_WARNING( "The field " << tag << " is ignored by the MSH writer, msh files do not support " << type << " data"); continue; } auto && data = mesh.getNodalData(tag); outfile << "$NodeData" << "\n"; outfile << "1" << "\n"; outfile << "\"" << tag << "\"\n"; outfile << "1\n0.0" << "\n"; outfile << "3\n0" << "\n"; outfile << data.getNbComponent() << "\n"; outfile << data.size() << "\n"; for (auto && d : enumerate(make_view(data, data.getNbComponent()))) { outfile << std::get<0>(d) + 1; for (auto && v : std::get<1>(d)) { outfile << " " << v; } outfile << "\n"; } outfile << "$EndNodeData" << "\n"; } } if (mesh.hasData(MeshDataType::_elemental)) { auto && tags = mesh.getTagNames(); for (auto && tag : tags) { auto && data = mesh.getElementalData(tag); if (mesh.getTypeCode(tag, MeshDataType::_elemental) != _tc_real or not data.isNodal()) continue; auto size = data.size(); if (size == 0) continue; auto && nb_components = data.getNbComponents(); auto nb_component = nb_components(*(data.elementTypes().begin())); outfile << "$ElementData" << "\n"; outfile << "1" << "\n"; outfile << "\"" << tag << "\"\n"; outfile << "1\n0.0" << "\n"; outfile << "3\n0" << "\n"; outfile << nb_component << "\n"; outfile << size << "\n"; Element element; for (auto type : data.elementTypes()) { element.type = type; for (auto && _ : enumerate(make_view(data(type), nb_components(type)))) { element.element = std::get<0>(_); outfile << element_to_msh_element[element]; for (auto && v : std::get<1>(_)) { outfile << " " << v; } outfile << "\n"; } } outfile << "$EndElementData" << "\n"; } } outfile.close(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh index 982793220..cd464b76d 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,468 +1,468 @@ /** * @file element_type_map.hh * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Tue Feb 20 2018 * * @brief storage class by element type * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_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(nb_component_functor); 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. */ template inline Stored & operator()(U && insertee, const SupportType & type, const GhostType & ghost_type = _not_ghost); public: /// 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; /// helper class to use in range for constructions 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)) {} + OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined)) {} ElementTypesIteratorHelper(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(ElementTypesIteratorHelper &&) = default; iterator begin(); iterator end(); 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; + ElementKind kind = _ek_not_defined) const; template ElementTypesIteratorHelper elementTypesImpl(const use_named_args_t & /*unused*/, pack &&... _pack) const; public: /*! * \param _spatial_dimension optional: filter for elements of given spatial * dimension * \param _ghost_type optional: filter for a certain ghost_type * \param _element_kind optional: filter for elements of given kind */ 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)...); } /*! 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*/ [[deprecated("Use elementTypes instead")]] 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 */ [[deprecated("Use elementTypes instead")]] inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; /*! 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; /// explicit copy void copy(const ElementTypeMapArray & other); /*! 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 */ 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); /** * get the size of the ElementTypeMapArray * @param[in] _spatial_dimension the dimension to consider (default: * _all_dimensions) * @param[in] _ghost_type (default: _not_ghost) * @param[in] _element_kind (default: _ek_not_defined) * @param[in] _all_ghost_types (default: false) **/ template UInt size(pack &&... _pack) const; bool isNodal() const { return is_nodal; } void isNodal(bool is_nodal) { this->is_nodal = is_nodal; } private: UInt sizeImpl(UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & kind) const; protected: /// name of the element type map: e.g. connectivity, grad_u ID name; /// Is the data stored by node of the element bool is_nodal{false}; }; /// 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 42e5e1b6b..776803f85 100644 --- a/src/mesh/element_type_map_tmpl.hh +++ b/src/mesh/element_type_map_tmpl.hh @@ -1,784 +1,784 @@ /** * @file element_type_map_tmpl.hh * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Tue Feb 20 2018 * * @brief implementation of template functions of the ElementTypeMap and * ElementTypeMapArray classes * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_static_if.hh" #include "element_type_map.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include "element_type_conversion.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #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 template inline Stored & ElementTypeMap:: operator()(U && insertee, 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::make_pair(type, std::forward(insertee))); 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 void ElementTypeMapArray::copy( const ElementTypeMapArray & other) { for (auto ghost_type : ghost_types) { for (auto type : this->elementTypes(_all_dimensions, ghost_types, _ek_not_defined)) { const auto & array_to_copy = other(type, ghost_type); auto & array = this->alloc(0, array_to_copy.getNbComponent(), type, ghost_type); array.copy(array_to_copy); } } } /* -------------------------------------------------------------------------- */ 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()) { auto id = this->id + ":" + aka::to_string(type) + ghost_id; tmp = &(Memory::alloc(id, size, nb_component, default_value)); 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 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 auto ElementTypeMap::ElementTypesIteratorHelper::begin() -> iterator { auto b = container.get().getData(ghost_type).begin(); auto e = container.get().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 iterator(b, e, dim, kind); } template auto ElementTypeMap::ElementTypesIteratorHelper::end() -> iterator { auto e = container.get().getData(ghost_type).end(); return iterator(e, e, dim, kind); } /* -------------------------------------------------------------------------- */ template auto ElementTypeMap::elementTypesImpl( UInt dim, GhostType ghost_type, ElementKind kind) const -> ElementTypesIteratorHelper { return ElementTypesIteratorHelper(*this, dim, ghost_type, kind); } /* -------------------------------------------------------------------------- */ template template auto ElementTypeMap::elementTypesImpl( const use_named_args_t & unused, pack &&... _pack) const -> ElementTypesIteratorHelper { return ElementTypesIteratorHelper(*this, unused, _pack...); } /* -------------------------------------------------------------------------- */ template inline auto ElementTypeMap::firstType( UInt dim, GhostType ghost_type, ElementKind kind) const -> type_iterator { return elementTypes(dim, ghost_type, kind).begin(); } /* -------------------------------------------------------------------------- */ template inline auto ElementTypeMap::lastType( UInt dim, GhostType ghost_type, ElementKind kind) const -> type_iterator { 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 ElementTypeMapArrayInitializer { protected: using CompFunc = std::function; public: ElementTypeMapArrayInitializer(const CompFunc & comp_func, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, - const ElementKind & element_kind = _ek_regular) + const ElementKind & element_kind = _ek_not_defined) : comp_func(comp_func), spatial_dimension(spatial_dimension), ghost_type(ghost_type), element_kind(element_kind) {} const GhostType & ghostType() const { return ghost_type; } virtual UInt nbComponent(const ElementType & type) const { return comp_func(type, ghostType()); } virtual bool isNodal() const { return false; } protected: CompFunc comp_func; UInt spatial_dimension; GhostType ghost_type; ElementKind element_kind; }; /* -------------------------------------------------------------------------- */ class MeshElementTypeMapArrayInitializer : public ElementTypeMapArrayInitializer { using CompFunc = ElementTypeMapArrayInitializer::CompFunc; public: MeshElementTypeMapArrayInitializer( const Mesh & mesh, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, - const ElementKind & element_kind = _ek_regular, + const ElementKind & element_kind = _ek_not_defined, bool with_nb_element = false, bool with_nb_nodes_per_element = false) : MeshElementTypeMapArrayInitializer( mesh, [nb_component](const ElementType &, const GhostType &) -> UInt { return nb_component; }, spatial_dimension, ghost_type, element_kind, with_nb_element, with_nb_nodes_per_element) {} MeshElementTypeMapArrayInitializer( const Mesh & mesh, const CompFunc & comp_func, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, - const ElementKind & element_kind = _ek_regular, + const ElementKind & element_kind = _ek_not_defined, bool with_nb_element = false, bool with_nb_nodes_per_element = false) : ElementTypeMapArrayInitializer(comp_func, spatial_dimension, 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 override { UInt res = ElementTypeMapArrayInitializer::nbComponent(type); if (with_nb_nodes_per_element) return (res * mesh.getNbNodesPerElement(type)); return res; } bool isNodal() const override { return with_nb_nodes_per_element; } protected: const Mesh & mesh; bool with_nb_element; bool with_nb_nodes_per_element; }; /* -------------------------------------------------------------------------- */ class FEEngineElementTypeMapArrayInitializer : public MeshElementTypeMapArrayInitializer { public: FEEngineElementTypeMapArrayInitializer( 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); + const ElementKind & element_kind = _ek_not_defined); FEEngineElementTypeMapArrayInitializer( const FEEngine & fe_engine, const ElementTypeMapArrayInitializer::CompFunc & nb_component, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, - const ElementKind & element_kind = _ek_regular); + const ElementKind & element_kind = _ek_not_defined); 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) { this->is_nodal = f.isNodal(); auto ghost_type = f.ghostType(); for (auto & type : f.elementTypes()) { if (not this->exists(type, ghost_type)) if (do_not_default) { auto & array = this->alloc(0, f.nbComponent(type), type, ghost_type); array.resize(f.size(type)); } else { this->alloc(f.size(type), f.nbComponent(type), type, ghost_type, default_value); } else { auto & array = this->operator()(type, ghost_type); if (not do_not_default) array.resize(f.size(type), default_value); else array.resize(f.size(type)); } } } /* -------------------------------------------------------------------------- */ /** * All parameters are named optionals * \param _nb_component a functor giving the number of components per * (ElementType, GhostType) pair or a scalar giving a unique number of * components * regardless of type * \param _spatial_dimension a filter for the elements of a specific dimension * \param _element_kind filter with element kind (_ek_regular, _ek_structural, * ...) * \param _with_nb_element allocate the arrays with the number of elements for * each * type in the mesh * \param _with_nb_nodes_per_element multiply the number of components by the * number of nodes per element * \param _default_value default inital value * \param _do_not_default do not initialize the allocated arrays * \param _ghost_type filter a type of ghost */ 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; auto functor = MeshElementTypeMapArrayInitializer( mesh, OPTIONAL_NAMED_ARG(nb_component, 1), OPTIONAL_NAMED_ARG(spatial_dimension, mesh.getSpatialDimension()), - ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_regular), + ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined), OPTIONAL_NAMED_ARG(with_nb_element, false), OPTIONAL_NAMED_ARG(with_nb_nodes_per_element, false)); this->initialize(functor, OPTIONAL_NAMED_ARG(default_value, T()), OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ /** * All parameters are named optionals * \param _nb_component a functor giving the number of components per * (ElementType, GhostType) pair or a scalar giving a unique number of * components * regardless of type * \param _spatial_dimension a filter for the elements of a specific dimension * \param _element_kind filter with element kind (_ek_regular, _ek_structural, * ...) * \param _default_value default inital value * \param _do_not_default do not initialize the allocated arrays * \param _ghost_type filter a specific ghost type * \param _all_ghost_types get all ghost types */ 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; auto functor = FEEngineElementTypeMapArrayInitializer( 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(element_kind, _ek_not_defined)); this->initialize(functor, 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); } /* -------------------------------------------------------------------------- */ template UInt ElementTypeMapArray::sizeImpl( UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & kind) const { UInt size = 0; for(auto && type : this->elementTypes(spatial_dimension, ghost_type, kind)) { size += this->operator()(type, ghost_type).size(); } return size; } /* -------------------------------------------------------------------------- */ template template UInt ElementTypeMapArray::size(pack &&... _pack) const { UInt size = 0; 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; size += sizeImpl(OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_not_defined)); } return size; } } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ */ diff --git a/src/mesh/group_manager.hh b/src/mesh/group_manager.hh index fb542763d..8879b03c0 100644 --- a/src/mesh/group_manager.hh +++ b/src/mesh/group_manager.hh @@ -1,356 +1,364 @@ /** * @file group_manager.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Feb 07 2018 * * @brief Stores information relevent to the notion of element and nodes * groups. * * @section LICENSE * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GROUP_MANAGER_HH__ #define __AKANTU_GROUP_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_iterators.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class ElementGroup; class NodeGroup; class Mesh; class Element; class ElementSynchronizer; template class CommunicationBufferTemplated; namespace dumper { class Field; } } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class GroupManager { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: #ifdef SWIGPYTHON public: using ElementGroups = std::map; using NodeGroups = std::map; private: #else using ElementGroups = std::map; using NodeGroups = std::map; #endif public: using GroupManagerTypeSet = std::set; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GroupManager(const Mesh & mesh, const ID & id = "group_manager", const MemoryID & memory_id = 0); virtual ~GroupManager(); /* ------------------------------------------------------------------------ */ /* Groups iterators */ /* ------------------------------------------------------------------------ */ public: using node_group_iterator = NodeGroups::iterator; using element_group_iterator = ElementGroups::iterator; using const_node_group_iterator = NodeGroups::const_iterator; using const_element_group_iterator = ElementGroups::const_iterator; #ifndef SWIG #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(group_type, function, \ param_in, param_out) \ inline BOOST_PP_CAT(BOOST_PP_CAT(const_, group_type), _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) const { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ }; \ \ inline BOOST_PP_CAT(group_type, _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ } #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(group_type, function) \ AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION( \ group_type, function, BOOST_PP_EMPTY(), BOOST_PP_EMPTY()) AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(element_group, find, const std::string & name, name); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(node_group, find, const std::string & name, name); #endif public: #ifndef SWIG decltype(auto) iterateNodeGroups() { return make_dereference_adaptor(make_values_adaptor(node_groups)); } decltype(auto) iterateNodeGroups() const { return make_dereference_adaptor(make_values_adaptor(node_groups)); } + + decltype(auto) iterateElementGroups() { + return make_dereference_adaptor(make_values_adaptor(element_groups)); + } + decltype(auto) iterateElementGroups() const { + return make_dereference_adaptor(make_values_adaptor(element_groups)); + } + #endif /* ------------------------------------------------------------------------ */ /* Clustering filter */ /* -------------------------------------------------------------------9+ ----- */ public: class ClusteringFilter { public: virtual bool operator()(const Element &) const { return true; } }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// create an empty node group NodeGroup & createNodeGroup(const std::string & group_name, bool replace_group = false); /// create an element group and the associated node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension = _all_dimensions, bool replace_group = false); /* ------------------------------------------------------------------------ */ /// renames an element group void renameElementGroup(const std::string & name, const std::string & new_name); /// renames a node group void renameNodeGroup(const std::string & name, const std::string & new_name); /// copy an existing element group void copyElementGroup(const std::string & name, const std::string & new_name); /// copy an existing node group void copyNodeGroup(const std::string & name, const std::string & new_name); /* ------------------------------------------------------------------------ */ /// create a node group from another node group but filtered template NodeGroup & createFilteredNodeGroup(const std::string & group_name, const NodeGroup & node_group, T & filter); /// destroy a node group void destroyNodeGroup(const std::string & group_name); /// create an element group from another element group but filtered template ElementGroup & createFilteredElementGroup(const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter); /// destroy an element group and the associated node group void destroyElementGroup(const std::string & group_name, bool destroy_node_group = false); /// destroy all element groups and the associated node groups void destroyAllElementGroups(bool destroy_node_groups = false); /// create a element group using an existing node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group); /// create groups based on values stored in a given mesh data template void createGroupsFromMeshData(const std::string & dataset_name); /// create boundaries group by a clustering algorithm \todo extend to parallel UInt createBoundaryGroupFromGeometry(); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, Mesh & mesh_facets, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); private: /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, const std::string & cluster_name_prefix, const ClusteringFilter & filter, Mesh & mesh_facets); public: /// Create an ElementGroup based on a NodeGroup void createElementGroupFromNodeGroup(const std::string & name, const std::string & node_group, UInt dimension = _all_dimensions); virtual void printself(std::ostream & stream, int indent = 0) const; /// this function insure that the group names are present on all processors /// /!\ it is a SMP call void synchronizeGroupNames(); /// register an elemental field to the given group name (overloading for /// ElementalPartionField) #ifndef SWIG template class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// ElementalField) template class ret_type, template class, bool> class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// MaterialInternalField) template class dump_type> dumper::Field * createElementalField(const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); template class ftype> dumper::Field * createNodalField(const ftype * field, const std::string & group_name, UInt padding_size = 0); template class ftype> dumper::Field * createStridedNodalField(const ftype * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size); protected: /// fill a buffer with all the group names void fillBufferWithGroupNames( CommunicationBufferTemplated & comm_buffer) const; /// take a buffer and create the missing groups localy void checkAndAddGroups(CommunicationBufferTemplated & buffer); /// register an elemental field to the given group name template inline dumper::Field * createElementalField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, const ElementTypeMap & nb_data_per_elem); /// register an elemental field to the given group name template inline dumper::Field * createElementalFilteredField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); #endif /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: - AKANTU_GET_MACRO(ElementGroups, element_groups, const ElementGroups &); + //AKANTU_GET_MACRO(ElementGroups, element_groups, const ElementGroups &); const ElementGroup & getElementGroup(const std::string & name) const; const NodeGroup & getNodeGroup(const std::string & name) const; ElementGroup & getElementGroup(const std::string & name); NodeGroup & getNodeGroup(const std::string & name); UInt getNbElementGroups(UInt dimension = _all_dimensions) const; UInt getNbNodeGroups() { return node_groups.size(); }; bool elementGroupExists(const std::string & name) { return element_groups.find(name) != element_groups.end(); } bool nodeGroupExists(const std::string & name) { return node_groups.find(name) != node_groups.end(); } private: template void renameGroup(GroupsType & groups, const std::string & name, const std::string & new_name); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id to create element and node groups ID id; /// memory_id to create element and node groups MemoryID memory_id; /// list of the node groups managed NodeGroups node_groups; /// list of the element groups managed ElementGroups element_groups; /// Mesh to which the element belongs const Mesh & mesh; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const GroupManager & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AKANTU_GROUP_MANAGER_HH__ */ diff --git a/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index 2eaaf3275..8f934046d 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,952 +1,952 @@ /** * @file heat_transfer_model.cc * * @author Guillaume Anciaux * @author Lucas Frerot * @author Emil Gallyamov * @author David Simon Kammer * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Tue Feb 20 2018 * * @brief Implementation of HeatTransferModel class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_model.hh" #include "dumpable_inline_impl.hh" #include "element_synchronizer.hh" #include "fe_engine_template.hh" #include "generalized_trapezoidal.hh" #include "group_manager_inline_impl.cc" #include "integrator_gauss.hh" #include "mesh.hh" #include "parser.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { namespace heat_transfer { namespace details { class ComputeRhoFunctor { public: ComputeRhoFunctor(const HeatTransferModel & model) : model(model){}; void operator()(Matrix & rho, const Element &) const { rho.set(model.getCapacity()*model.getDensity()); } 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) { AKANTU_DEBUG_IN(); conductivity = Matrix(this->spatial_dimension, this->spatial_dimension); this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, _gst_htm_temperature); this->registerSynchronizer(synchronizer, _gst_htm_gradient_temperature); } registerFEEngineObject(id + ":fem", mesh, spatial_dimension); #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("heat_transfer", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif this->registerParam("conductivity", conductivity, _pat_parsmod); this->registerParam("conductivity_variation", conductivity_variation, 0., _pat_parsmod); this->registerParam("temperature_reference", T_ref, 0., _pat_parsmod); this->registerParam("capacity", capacity, _pat_parsmod); this->registerParam("density", density, _pat_parsmod); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initModel() { auto & fem = this->getFEEngine(); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); temperature_on_qpoints.initialize(fem, _nb_component = 1); temperature_gradient.initialize(fem, _nb_component = spatial_dimension); conductivity_on_qpoints.initialize( fem, _nb_component = spatial_dimension * spatial_dimension); k_gradt_on_qpoints.initialize(fem, _nb_component = spatial_dimension); } /* -------------------------------------------------------------------------- */ FEEngine & HeatTransferModel::getFEEngineBoundary(const ID & name) { return dynamic_cast(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(); heat_transfer::details::ComputeRhoFunctor compute_rho(*this); - for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { + for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldLumped(compute_rho, "M", "temperature", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType HeatTransferModel::getMatrixType(const ID & matrix_id) { if (matrix_id == "K" or matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleConductivityMatrix(); } else if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacity(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacityLumped(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleResidual() { AKANTU_DEBUG_IN(); this->assembleInternalHeatRate(); this->getDOFManager().assembleToResidual("temperature", *this->external_heat_rate, 1); this->getDOFManager().assembleToResidual("temperature", *this->internal_heat_rate, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::predictor() { ++temperature_release; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped() { AKANTU_DEBUG_IN(); if (!this->getDOFManager().hasLumpedMatrix("M")) { this->getDOFManager().getNewLumpedMatrix("M"); } this->getDOFManager().clearLumpedMatrix("M"); assembleCapacityLumped(_not_ghost); assembleCapacityLumped(_ghost); need_to_reassemble_capacity_lumped = false; 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; } else { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["temperature"] = _ist_backward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleConductivityMatrix() { AKANTU_DEBUG_IN(); this->computeConductivityOnQuadPoints(_not_ghost); if (conductivity_release[_not_ghost] == conductivity_matrix_release) return; if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().clearMatrix("K"); switch (mesh.getSpatialDimension()) { case 1: this->assembleConductivityMatrix<1>(_not_ghost); break; case 2: this->assembleConductivityMatrix<2>(_not_ghost); break; case 3: this->assembleConductivityMatrix<3>(_not_ghost); break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::assembleConductivityMatrix( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine(); - for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type)) { + for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto nb_element = mesh.getNbElement(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); fem.computeBtDB(conductivity_on_qpoints(type, 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); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "temperature", *K_e, type, ghost_type, _symmetric); } conductivity_matrix_release = conductivity_release[ghost_type]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeConductivityOnQuadPoints( const GhostType & ghost_type) { // if already computed once check if need to compute if (not initial_conductivity[ghost_type]) { // if temperature did not change, condictivity will not vary if (temperature_release == conductivity_release[ghost_type]) return; // if conductivity_variation is 0 no need to recompute if (conductivity_variation == 0.) return; } - for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { + for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & temperature_interpolated = temperature_on_qpoints(type, ghost_type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type, ghost_type); auto & cond = conductivity_on_qpoints(type, ghost_type); for (auto && tuple : zip(make_view(cond, spatial_dimension, spatial_dimension), temperature_interpolated)) { auto & C = std::get<0>(tuple); auto & T = std::get<1>(tuple); C = conductivity; Matrix variation(spatial_dimension, spatial_dimension, conductivity_variation * (T - T_ref)); // @TODO: Guillaume are you sure ? why due you compute variation then ? C += conductivity_variation; } } conductivity_release[ghost_type] = temperature_release; initial_conductivity[ghost_type] = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeKgradT(const GhostType & ghost_type) { computeConductivityOnQuadPoints(ghost_type); - for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) { + for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & gradient = temperature_gradient(type, ghost_type); this->getFEEngine().gradientOnIntegrationPoints(*temperature, gradient, 1, type, ghost_type); for (auto && values : zip(make_view(conductivity_on_qpoints(type, ghost_type), spatial_dimension, spatial_dimension), make_view(gradient, spatial_dimension), make_view(k_gradt_on_qpoints(type, ghost_type), spatial_dimension))) { const auto & C = std::get<0>(values); const auto & BT = std::get<1>(values); auto & k_BT = std::get<2>(values); k_BT.mul(C, BT); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleInternalHeatRate() { AKANTU_DEBUG_IN(); this->internal_heat_rate->clear(); this->synchronize(_gst_htm_temperature); auto & fem = this->getFEEngine(); for (auto ghost_type : ghost_types) { // compute k \grad T computeKgradT(ghost_type); - for (auto type : mesh.elementTypes(spatial_dimension, ghost_type)) { + for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto & k_gradt_on_qpoints_vect = k_gradt_on_qpoints(type, ghost_type); UInt nb_quad_points = k_gradt_on_qpoints_vect.size(); Array bt_k_gT(nb_quad_points, nb_nodes_per_element); fem.computeBtD(k_gradt_on_qpoints_vect, bt_k_gT, type, ghost_type); UInt nb_elements = mesh.getNbElement(type, ghost_type); Array int_bt_k_gT(nb_elements, nb_nodes_per_element); fem.integrate(bt_k_gT, int_bt_k_gT, nb_nodes_per_element, type, ghost_type); this->getDOFManager().assembleElementalArrayLocalArray( int_bt_k_gT, *this->internal_heat_rate, type, ghost_type, -1); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real el_size; Real min_el_size = std::numeric_limits::max(); Real conductivitymax = conductivity(0, 0); // get the biggest parameter from k11 until k33// for (UInt i = 0; i < spatial_dimension; i++) for (UInt j = 0; j < spatial_dimension; j++) conductivitymax = std::max(conductivity(i, j), conductivitymax); - for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost)) { + for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array coord(0, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, _not_ghost); auto el_coord = coord.begin(spatial_dimension, nb_nodes_per_element); UInt nb_element = mesh.getNbElement(type); for (UInt el = 0; el < nb_element; ++el, ++el_coord) { el_size = getFEEngine().getElementInradius(*el_coord, type); min_el_size = std::min(min_el_size, el_size); } AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size << " and the max conductivity is : " << conductivitymax); } Real min_dt = 2. * min_el_size * min_el_size / 4. * density * capacity / conductivitymax; mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("heat_transfer").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void HeatTransferModel::readMaterials() { auto sect = this->getParserSection(); if (not std::get<1>(sect)) { const auto & section = std::get<0>(sect); this->parseSection(section); } conductivity_on_qpoints.set(conductivity); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); readMaterials(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacity() { AKANTU_DEBUG_IN(); auto ghost_type = _not_ghost; this->getDOFManager().clearMatrix("M"); auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor rho_functor(*this); - for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type)) { + for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldMatrix(rho_functor, "M", "temperature", this->getDOFManager(), type, ghost_type); } need_to_reassemble_capacity = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeRho(Array & rho, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); FEEngine & fem = this->getFEEngine(); UInt nb_element = mesh.getNbElement(type, ghost_type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); rho.resize(nb_element * nb_quadrature_points); rho.set(this->capacity); // Real * rho_1_val = rho.storage(); // /// compute @f$ rho @f$ for each nodes of each element // for (UInt el = 0; el < nb_element; ++el) { // for (UInt n = 0; n < nb_quadrature_points; ++n) { // *rho_1_val++ = this->capacity; // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::computeThermalEnergyByNode() { AKANTU_DEBUG_IN(); Real ethermal = 0.; for (auto && pair : enumerate(make_view( *internal_heat_rate, internal_heat_rate->getNbComponent()))) { auto n = std::get<0>(pair); auto & heat_rate = std::get<1>(pair); Real heat = 0.; bool is_local_node = mesh.isLocalOrMasterNode(n); bool count_node = is_local_node; for (UInt i = 0; i < heat_rate.size(); ++i) { if (count_node) heat += heat_rate[i] * time_step; } ethermal += heat; } mesh.getCommunicator().allReduce(ethermal, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ethermal; } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::getThermalEnergy( iterator Eth, Array::const_iterator T_it, 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; auto & fem = getFEEngine(); - for (auto && type : mesh.elementTypes(spatial_dimension)) { + for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_element = mesh.getNbElement(type, _not_ghost); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, _not_ghost); Array Eth_per_quad(nb_element * nb_quadrature_points, 1); auto & temperature_interpolated = temperature_on_qpoints(type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type); auto T_it = temperature_interpolated.begin(); auto T_end = temperature_interpolated.end(); getThermalEnergy(Eth_per_quad.begin(), T_it, T_end); Eth += fem.integrate(Eth_per_quad, type); } return Eth; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id) { AKANTU_DEBUG_IN(); Real energy = 0; if (id == "thermal") energy = getThermalEnergy(); // reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id, 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) { if (field_name == "capacity_lumped"){ AKANTU_EXCEPTION("Capacity lumped is a nodal field now stored in the DOF manager." "Therefore it cannot be used by a dumper anymore"); } std::map *> real_nodal_fields; real_nodal_fields["temperature"] = temperature; 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["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 /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void HeatTransferModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbData(const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = indexes.size(); switch (tag) { case _gst_htm_temperature: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case _gst_htm_temperature: { buffer << (*temperature)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case _gst_htm_temperature: { buffer >> (*temperature)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; Array::const_iterator it = elements.begin(); Array::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case _gst_htm_temperature: { size += nb_nodes_per_element * sizeof(Real); // temperature break; } case _gst_htm_gradient_temperature: { // temperature gradient size += getNbIntegrationPoints(elements) * spatial_dimension * sizeof(Real); size += nb_nodes_per_element * sizeof(Real); // nodal temperatures break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { switch (tag) { case _gst_htm_temperature: { packNodalDataHelper(*temperature, buffer, elements, mesh); break; } case _gst_htm_gradient_temperature: { packElementalDataHelper(temperature_gradient, buffer, elements, true, getFEEngine()); packNodalDataHelper(*temperature, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { switch (tag) { case _gst_htm_temperature: { unpackNodalDataHelper(*temperature, buffer, elements, mesh); break; } case _gst_htm_gradient_temperature: { unpackElementalDataHelper(temperature_gradient, buffer, elements, true, getFEEngine()); unpackNodalDataHelper(*temperature, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ } // akantu diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index bdaef7e82..584eed030 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1365 +1,1365 @@ /** * @file material.cc * * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_material, id), is_init(false), fem(model.getFEEngine()), finite_deformation(false), name(""), model(model), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id, this->memory_id), stress("stress", *this), eigengradu("eigen_grad_u", *this), gradu("grad_u", *this), green_strain("green_strain", *this), piola_kirchhoff_2("piola_kirchhoff_2", *this), potential_energy("potential_energy", *this), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse coordinates", *this), interpolation_points_matrices("interpolation points matrices", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(model.getMesh(), - _spatial_dimension = spatial_dimension); + _spatial_dimension = spatial_dimension, + _element_kind = _ek_regular); // model.getMesh().initElementTypeMapArray(element_filter, 1, // spatial_dimension, // false, _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_material, id), is_init(false), fem(fe_engine), finite_deformation(false), name(""), model(model), spatial_dimension(dim), element_filter("element_filter", id, this->memory_id), stress("stress", *this, dim, fe_engine, this->element_filter), eigengradu("eigen_grad_u", *this, dim, fe_engine, this->element_filter), gradu("gradu", *this, dim, fe_engine, this->element_filter), green_strain("green_strain", *this, dim, fe_engine, this->element_filter), piola_kirchhoff_2("piola_kirchhoff_2", *this, dim, fe_engine, this->element_filter), potential_energy("potential_energy", *this, dim, fe_engine, this->element_filter), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse_coordinates", *this, dim, fe_engine, this->element_filter), interpolation_points_matrices("interpolation points matrices", *this, dim, fe_engine, this->element_filter) { AKANTU_DEBUG_IN(); - element_filter.initialize(mesh, _spatial_dimension = spatial_dimension); + element_filter.initialize(mesh, _spatial_dimension = spatial_dimension, + _element_kind = _ek_regular); // mesh.initElementTypeMapArray(element_filter, 1, spatial_dimension, false, // _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() = default; /* -------------------------------------------------------------------------- */ void Material::initialize() { registerParam("rho", rho, Real(0.), _pat_parsable | _pat_modifiable, "Density"); registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("finite_deformation", finite_deformation, false, _pat_parsable | _pat_readable, "Is finite deformation"); registerParam("inelastic_deformation", inelastic_deformation, false, _pat_internal, "Is inelastic deformation"); /// allocate gradu stress for local elements eigengradu.initialize(spatial_dimension * spatial_dimension); gradu.initialize(spatial_dimension * spatial_dimension); stress.initialize(spatial_dimension * spatial_dimension); potential_energy.initialize(1); this->model.registerEventHandler(*this); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); if (finite_deformation) { this->piola_kirchhoff_2.initialize(spatial_dimension * spatial_dimension); if (use_previous_stress) this->piola_kirchhoff_2.initializeHistory(); this->green_strain.initialize(spatial_dimension * spatial_dimension); } if (use_previous_stress) this->stress.initializeHistory(); if (use_previous_gradu) this->gradu.initializeHistory(); this->resizeInternals(); is_init = true; updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::savePreviousState() { AKANTU_DEBUG_IN(); for (auto pair : internal_vectors_real) if (pair.second->hasHistory()) pair.second->saveCurrentValues(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::restorePreviousState() { AKANTU_DEBUG_IN(); for (auto pair : internal_vectors_real) if (pair.second->hasHistory()) pair.second->restorePreviousValues(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ // void Material::updateResidual(GhostType ghost_type) { // AKANTU_DEBUG_IN(); // computeAllStresses(ghost_type); // assembleResidual(ghost_type); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); if (!finite_deformation) { auto & internal_force = const_cast &>(model.getInternalForce()); // Mesh & mesh = fem.getMesh(); for (auto && type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) continue; const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by /// @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Array * sigma_dphi_dx = new Array(nb_element * nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); fem.computeBtD(stress(type, ghost_type), *sigma_dphi_dx, type, ghost_type, elem_filter); /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by * @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ Array * int_sigma_dphi_dx = new Array(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); fem.integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, type, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble model.getDOFManager().assembleElementalArrayLocalArray( *int_sigma_dphi_dx, internal_force, type, ghost_type, -1, elem_filter); delete int_sigma_dphi_dx; } } else { switch (spatial_dimension) { case 1: this->assembleInternalForces<1>(ghost_type); break; case 2: this->assembleInternalForces<2>(ghost_type); break; case 3: this->assembleInternalForces<3>(ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the gradu * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (const auto & type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, spatial_dimension, type, ghost_type, elem_filter); gradu_vect -= eigengradu(type, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllCauchyStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(finite_deformation, "The Cauchy stress can only be " "computed if you are working in " "finite deformation."); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { switch (spatial_dimension) { case 1: this->computeCauchyStress<1>(type, ghost_type); break; case 2: this->computeCauchyStress<2>(type, ghost_type); break; case 3: this->computeCauchyStress<3>(type, ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeCauchyStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::matrix_iterator gradu_it = this->gradu(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator gradu_end = this->gradu(el_type, ghost_type).end(dim, dim); Array::matrix_iterator piola_it = this->piola_kirchhoff_2(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator stress_it = this->stress(el_type, ghost_type).begin(dim, dim); Matrix F_tensor(dim, dim); for (; gradu_it != gradu_end; ++gradu_it, ++piola_it, ++stress_it) { Matrix & grad_u = *gradu_it; Matrix & piola = *piola_it; Matrix & sigma = *stress_it; gradUToF(grad_u, F_tensor); this->computeCauchyStressOnQuad(F_tensor, piola, sigma); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & displacement = model.getDisplacement(); // resizeInternalArray(gradu); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(displacement, gradu_vect, spatial_dimension, type, ghost_type, elem_filter); setToSteadyState(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { if (finite_deformation) { switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL<1>(type, ghost_type); assembleStiffnessMatrixL2<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrixNL<2>(type, ghost_type); assembleStiffnessMatrixL2<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrixNL<3>(type, ghost_type); assembleStiffnessMatrixL2<3>(type, ghost_type); break; } } } else { switch (spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(type, ghost_type); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) { AKANTU_DEBUG_OUT(); return; } // const Array & shapes_derivatives = // fem.getShapesDerivatives(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); fem.computeBtDB(*tangent_stiffness_matrix, *bt_d_b, 4, type, ghost_type, elem_filter); delete tangent_stiffness_matrix; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); // Array & gradu_vect = delta_gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); fem.filterElementalData(fem.getMesh(), shapes_derivatives, *shapes_derivatives_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_s_b_size = dim * nb_nodes_per_element; Array * bt_s_b = new Array(nb_element * nb_quadrature_points, bt_s_b_size * bt_s_b_size, "B^t*D*B"); UInt piola_matrix_size = getCauchyStressMatrixSize(dim); Matrix B(piola_matrix_size, bt_s_b_size); Matrix Bt_S(bt_s_b_size, piola_matrix_size); Matrix S(piola_matrix_size, piola_matrix_size); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); auto Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); auto piola_it = piola_kirchhoff_2(type, ghost_type).begin(dim, dim); for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) { auto & Bt_S_B = *Bt_S_B_it; const auto & Piola_kirchhoff_matrix = *piola_it; setCauchyStressMatrix(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.template mul(B, S); Bt_S_B.template mul(Bt_S, B); } delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_s_b_size * bt_s_b_size, "K_e"); fem.integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, elem_filter); delete bt_s_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); fem.filterElementalData(fem.getMesh(), shapes_derivatives, *shapes_derivatives_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); auto shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); auto Bt_D_B_it = bt_d_b->begin(bt_d_b_size, bt_d_b_size); auto grad_u_it = gradu_vect.begin(dim, dim); auto D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); auto D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) { const auto & grad_u = *grad_u_it; const auto & D = *D_it; auto & Bt_D_B = *Bt_D_B_it; // transferBMatrixToBL1 (*shapes_derivatives_filtered_it, B, // nb_nodes_per_element); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element); B += B2; Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); Mesh & mesh = fem.getMesh(); - for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { + for (auto type : element_filter.elementTypes(_ghost_type = ghost_type)) { const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array * shapesd_filtered = new Array( nb_element, size_of_shapes_derivatives, "filtered shapesd"); fem.filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); // Set stress vectors UInt stress_size = getTangentStiffnessVoigtSize(dim); // Set matrices B and BNL* UInt bt_s_size = dim * nb_nodes_per_element; auto * bt_s = new Array(nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); auto grad_u_it = this->gradu(type, ghost_type).begin(dim, dim); auto grad_u_end = this->gradu(type, ghost_type).end(dim, dim); auto stress_it = this->piola_kirchhoff_2(type, ghost_type).begin(dim, dim); shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1); Matrix S_vect(stress_size, 1); Matrix B_tensor(stress_size, bt_s_size); Matrix B2_tensor(stress_size, bt_s_size); for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) { auto & grad_u = *grad_u_it; auto & r_it = *bt_s_it; auto & S_it = *stress_it; setCauchyStressArray(S_it, S_vect); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element); B_tensor += B2_tensor; r_it.template mul(B_tensor, S_vect); } delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * r_e = new Array(nb_element, bt_s_size, "r_e"); fem.integrate(*bt_s, *r_e, bt_s_size, type, ghost_type, elem_filter); delete bt_s; model.getDOFManager().assembleElementalArrayLocalArray( *r_e, internal_force, type, ghost_type, -1., elem_filter); delete r_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElements() { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { computePotentialEnergy(type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType, GhostType) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElements(); /// integrate the potential energy for each type of elements for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { epot += fem.integrate(potential_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy(ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real epot = 0.; Vector epot_on_quad_points(fem.getNbIntegrationPoints(type)); computePotentialEnergyByElement(type, index, epot_on_quad_points); epot = fem.integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(const std::string & type) { AKANTU_DEBUG_IN(); if (type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(const std::string & energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "potential") return getPotentialEnergy(type, index); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation( const ElementTypeMapArray & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); this->fem.initElementalFieldInterpolationFromIntegrationPoints( interpolation_points_coordinates, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, &(this->element_filter)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type) { this->fem.interpolateElementalFieldFromIntegrationPoints( this->stress, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, result, ghost_type, &(this->element_filter)); } /* -------------------------------------------------------------------------- */ void Material::interpolateStressOnFacets( ElementTypeMapArray & result, ElementTypeMapArray & by_elem_result, const GhostType ghost_type) { interpolateStress(by_elem_result, ghost_type); UInt stress_size = this->stress.getNbComponent(); const Mesh & mesh = this->model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_fil = element_filter(type, ghost_type); Array & by_elem_res = by_elem_result(type, ghost_type); UInt nb_element = elem_fil.size(); UInt nb_element_full = this->model.getMesh().getNbElement(type, ghost_type); UInt nb_interpolation_points_per_elem = by_elem_res.size() / nb_element_full; const Array & facet_to_element = mesh_facets.getSubelementToElement(type, ghost_type); ElementType type_facet = Mesh::getFacetType(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); UInt nb_quad_per_facet = nb_interpolation_points_per_elem / nb_facet_per_elem; Element element_for_comparison{type, 0, ghost_type}; const Array> * element_to_facet = nullptr; GhostType current_ghost_type = _casper; Array * result_vec = nullptr; Array::const_matrix_iterator result_it = by_elem_res.begin_reinterpret( stress_size, nb_interpolation_points_per_elem, nb_element_full); for (UInt el = 0; el < nb_element; ++el) { UInt global_el = elem_fil(el); element_for_comparison.element = global_el; for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element facet_elem = facet_to_element(global_el, f); UInt global_facet = facet_elem.element; if (facet_elem.ghost_type != current_ghost_type) { current_ghost_type = facet_elem.ghost_type; element_to_facet = &mesh_facets.getElementToSubelement( type_facet, current_ghost_type); result_vec = &result(type_facet, current_ghost_type); } bool is_second_element = (*element_to_facet)(global_facet)[0] != element_for_comparison; for (UInt q = 0; q < nb_quad_per_facet; ++q) { Vector result_local(result_vec->storage() + (global_facet * nb_quad_per_facet + q) * result_vec->getNbComponent() + is_second_element * stress_size, stress_size); const Matrix & result_tmp(result_it[global_el]); result_local = result_tmp(f * nb_quad_per_facet + q); } } } } } /* -------------------------------------------------------------------------- */ template const Array & Material::getArray(const ID & /*vect_id*/, const ElementType & /*type*/, const GhostType & /*ghost_type*/) const { AKANTU_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template Array & Material::getArray(const ID & /*vect_id*/, const ElementType & /*type*/, const GhostType & /*ghost_type*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << " (" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << " (" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << " (" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template -const InternalField & Material::getInternal([[gnu::unused]] - const ID & int_id) const { +const InternalField & +Material::getInternal([[gnu::unused]] const ID & int_id) const { AKANTU_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template -InternalField & Material::getInternal([[gnu::unused]] - const ID & int_id) { +InternalField & Material::getInternal([[gnu::unused]] const ID & int_id) { AKANTU_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ void Material::addElements(const Array & elements_to_add) { AKANTU_DEBUG_IN(); UInt mat_id = model.getInternalIndexFromID(getID()); Array::const_iterator el_begin = elements_to_add.begin(); Array::const_iterator el_end = elements_to_add.end(); for (; el_begin != el_end; ++el_begin) { const Element & element = *el_begin; Array & mat_indexes = model.getMaterialByElement(element.type, element.ghost_type); Array & mat_loc_num = model.getMaterialLocalNumbering(element.type, element.ghost_type); UInt index = this->addElement(element.type, element.element, element.ghost_type); mat_indexes(element.element) = mat_id; mat_loc_num(element.element) = index; } this->resizeInternals(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::removeElements(const Array & elements_to_remove) { AKANTU_DEBUG_IN(); Array::const_iterator el_begin = elements_to_remove.begin(); Array::const_iterator el_end = elements_to_remove.end(); if (el_begin == el_end) return; ElementTypeMapArray material_local_new_numbering( "remove mat filter elem", getID(), getMemoryID()); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (auto & type : element_filter.elementTypes(_ghost_type = ghost_type, _element_kind = _ek_not_defined)) { element.type = type; Array & elem_filter = this->element_filter(type, ghost_type); Array & mat_loc_num = this->model.getMaterialLocalNumbering(type, ghost_type); if (!material_local_new_numbering.exists(type, ghost_type)) material_local_new_numbering.alloc(elem_filter.size(), 1, type, ghost_type); Array & mat_renumbering = material_local_new_numbering(type, ghost_type); UInt nb_element = elem_filter.size(); Array elem_filter_tmp; UInt new_id = 0; for (UInt el = 0; el < nb_element; ++el) { element.element = elem_filter(el); if (std::find(el_begin, el_end, element) == el_end) { elem_filter_tmp.push_back(element.element); mat_renumbering(el) = new_id; mat_loc_num(element.element) = new_id; ++new_id; } else { mat_renumbering(el) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.size()); elem_filter.copy(elem_filter_tmp); } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resizeInternals() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::onElementsAdded(const Array &, const NewElementsEvent &) { this->resizeInternals(); } /* -------------------------------------------------------------------------- */ void Material::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, [[gnu::unused]] const RemovedElementsEvent & event) { UInt my_num = model.getInternalIndexFromID(getID()); ElementTypeMapArray material_local_new_numbering( "remove mat filter elem", getID(), getMemoryID()); auto el_begin = element_list.begin(); auto el_end = element_list.end(); for (auto && gt : ghost_types) { for (auto && type : new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) { if (not element_filter.exists(type, gt) || element_filter(type, gt).size() == 0) continue; auto & elem_filter = element_filter(type, gt); auto & mat_indexes = this->model.getMaterialByElement(type, gt); auto & mat_loc_num = this->model.getMaterialLocalNumbering(type, gt); auto nb_element = this->model.getMesh().getNbElement(type, gt); // all materials will resize of the same size... mat_indexes.resize(nb_element); mat_loc_num.resize(nb_element); if (!material_local_new_numbering.exists(type, gt)) material_local_new_numbering.alloc(elem_filter.size(), 1, type, gt); auto & mat_renumbering = material_local_new_numbering(type, gt); const auto & renumbering = new_numbering(type, gt); Array elem_filter_tmp; UInt ni = 0; Element el{type, 0, gt}; for (UInt i = 0; i < elem_filter.size(); ++i) { el.element = elem_filter(i); if (std::find(el_begin, el_end, el) == el_end) { UInt new_el = renumbering(el.element); AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A not removed element as been badly renumbered"); elem_filter_tmp.push_back(new_el); mat_renumbering(i) = ni; mat_indexes(new_el) = my_num; mat_loc_num(new_el) = ni; ++ni; } else { mat_renumbering(i) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.size()); elem_filter.copy(elem_filter_tmp); } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); } /* -------------------------------------------------------------------------- */ void Material::beforeSolveStep() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::afterSolveStep() { for (auto & type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { this->updateEnergies(type, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDamageIteration() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onDamageUpdate() { for (auto & type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { this->updateEnergiesAfterDamage(type, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDump() { if (this->isFiniteDeformation()) this->computeAllCauchyStresses(_not_ghost); } /* -------------------------------------------------------------------------- */ void Material::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /// extrapolate internal values void Material::extrapolateInternal(const ID & id, const Element & element, - [[gnu::unused]] - const Matrix & point, + [[gnu::unused]] const Matrix & point, Matrix & extrapolated) { if (this->isInternal(id, element.kind())) { UInt nb_element = this->element_filter(element.type, element.ghost_type).size(); const ID name = this->getID() + ":" + id; UInt nb_quads = this->internal_vectors_real[name]->getFEEngine().getNbIntegrationPoints( element.type, element.ghost_type); const Array & internal = this->getArray(id, element.type, element.ghost_type); UInt nb_component = internal.getNbComponent(); Array::const_matrix_iterator internal_it = internal.begin_reinterpret(nb_component, nb_quads, nb_element); Element local_element = this->convertToLocalElement(element); /// instead of really extrapolating, here the value of the first GP /// is copied into the result vector. This works only for linear /// elements /// @todo extrapolate!!!! AKANTU_DEBUG_WARNING("This is a fix, values are not truly extrapolated"); const Matrix & values = internal_it[local_element.element]; UInt index = 0; Vector tmp(nb_component); for (UInt j = 0; j < values.cols(); ++j) { tmp = values(j); if (tmp.norm() > 0) { index = j; break; } } for (UInt i = 0; i < extrapolated.size(); ++i) { extrapolated(i) = values(index); } } else { Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); extrapolated = default_values; } } /* -------------------------------------------------------------------------- */ void Material::applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const GhostType ghost_type) { for (auto && type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { if (!element_filter(type, ghost_type).size()) continue; auto eigen_it = this->eigengradu(type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto eigen_end = this->eigengradu(type, ghost_type) .end(spatial_dimension, spatial_dimension); for (; eigen_it != eigen_end; ++eigen_it) { auto & current_eigengradu = *eigen_it; current_eigengradu = prescribed_eigen_grad_u; } } } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_mass.cc b/src/model/solid_mechanics/solid_mechanics_model_mass.cc index 19cdca16b..3b5eeb42c 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_mass.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_mass.cc @@ -1,150 +1,150 @@ /** * @file solid_mechanics_model_mass.cc * * @author Nicolas Richart * * @date creation: Tue Oct 05 2010 * @date last modification: Wed Nov 08 2017 * * @brief function handling mass computation * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "integrator_gauss.hh" #include "material.hh" #include "model_solver.hh" #include "shape_lagrange.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class ComputeRhoFunctor { public: explicit ComputeRhoFunctor(const SolidMechanicsModel & model) : model(model){}; void operator()(Matrix & rho, const Element & element) const { const Array & mat_indexes = model.getMaterialByElement(element.type, element.ghost_type); Real mat_rho = model.getMaterial(mat_indexes(element.element)).getParam("rho"); rho.set(mat_rho); } private: const SolidMechanicsModel & model; }; /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped() { AKANTU_DEBUG_IN(); if (not need_to_reassemble_lumped_mass) return; this->allocNodalField(this->mass, spatial_dimension, "mass"); mass->clear(); if (!this->getDOFManager().hasLumpedMatrix("M")) { this->getDOFManager().getNewLumpedMatrix("M"); } this->getDOFManager().clearLumpedMatrix("M"); assembleMassLumped(_not_ghost); assembleMassLumped(_ghost); this->getDOFManager().getLumpedMatrixPerDOFs("displacement", "M", *(this->mass)); /// for not connected nodes put mass to one in order to avoid #if !defined(AKANTU_NDEBUG) bool has_unconnected_nodes = false; auto mass_it = mass->begin_reinterpret(mass->size() * mass->getNbComponent()); auto mass_end = mass->end_reinterpret(mass->size() * mass->getNbComponent()); for (; mass_it != mass_end; ++mass_it) { if (std::abs(*mass_it) < std::numeric_limits::epsilon() || Math::isnan(*mass_it)) { has_unconnected_nodes = true; break; } } if (has_unconnected_nodes) AKANTU_DEBUG_WARNING("There are nodes that seem to not be connected to any " "elements, beware that they have lumped mass of 0."); #endif this->synchronize(_gst_smm_mass); need_to_reassemble_lumped_mass = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass() { AKANTU_DEBUG_IN(); if (not need_to_reassemble_mass) return; this->getDOFManager().clearMatrix("M"); assembleMass(_not_ghost); need_to_reassemble_mass = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem = getFEEngineClass(); ComputeRhoFunctor compute_rho(*this); - for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { + for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldLumped(compute_rho, "M", "displacement", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem = getFEEngineClass(); ComputeRhoFunctor compute_rho(*this); - for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { + for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldMatrix(compute_rho, "M", "displacement", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/test/test_model/patch_tests/patch_test_linear_fixture.hh b/test/test_model/patch_tests/patch_test_linear_fixture.hh index e866bf5e7..aa1f50b71 100644 --- a/test/test_model/patch_tests/patch_test_linear_fixture.hh +++ b/test/test_model/patch_tests/patch_test_linear_fixture.hh @@ -1,183 +1,183 @@ /** * @file patch_test_linear_fixture.hh * * @author Nicolas Richart * * @date creation: Tue Jan 30 2018 * @date last modification: Wed Jan 31 2018 * * @brief Fixture for linear patch tests * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "mesh_utils.hh" #include "model.hh" #include "test_gtest_utils.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PATCH_TEST_LINEAR_FIXTURE_HH__ #define __AKANTU_PATCH_TEST_LINEAR_FIXTURE_HH__ //#define DEBUG_TEST using namespace akantu; template class TestPatchTestLinear : public ::testing::Test { public: static constexpr ElementType type = type_::value; static constexpr size_t dim = ElementClass::getSpatialDimension(); virtual void SetUp() { mesh = std::make_unique(dim); mesh->read(aka::to_string(type) + ".msh"); MeshUtils::buildFacets(*mesh); mesh->createBoundaryGroupFromGeometry(); model = std::make_unique(*mesh); } virtual void TearDown() { model.reset(nullptr); mesh.reset(nullptr); } virtual void initModel(const AnalysisMethod & method, const std::string & material_file) { debug::setDebugLevel(dblError); getStaticParser().parse(material_file); this->model->initFull(_analysis_method = method); this->applyBC(); if (method != _static) this->model->setTimeStep(0.8 * this->model->getStableTimeStep()); } virtual void applyBC() { auto & boundary = this->model->getBlockedDOFs(); - for (auto & eg : mesh->getElementGroups()) { - for (const auto & node : eg.second->getNodeGroup()) { + for (auto & eg : mesh->iterateElementGroups()) { + for (const auto & node : eg.getNodeGroup()) { for (UInt s = 0; s < boundary.getNbComponent(); ++s) { boundary(node, s) = true; } } } } virtual void applyBConDOFs(const Array & dofs) { const auto & coordinates = this->mesh->getNodes(); - for (auto & eg : this->mesh->getElementGroups()) { - for (const auto & node : eg.second->getNodeGroup()) { + for (auto & eg : this->mesh->iterateElementGroups()) { + for (const auto & node : eg.getNodeGroup()) { this->setLinearDOF(dofs.begin(dofs.getNbComponent())[node], coordinates.begin(this->dim)[node]); } } } template Matrix prescribed_gradient(const V & dof) { Matrix gradient(dof.getNbComponent(), dim); for (UInt i = 0; i < gradient.rows(); ++i) { for (UInt j = 0; j < gradient.cols(); ++j) { gradient(i, j) = alpha(i, j + 1); } } return gradient; } template void checkGradient(const Gradient & gradient, const DOFs & dofs) { auto pgrad = prescribed_gradient(dofs); for (auto & grad : make_view(gradient, gradient.getNbComponent() / dim, dim)) { auto diff = grad - pgrad; auto gradient_error = diff.template norm() / grad.template norm(); EXPECT_NEAR(0, gradient_error, gradient_tolerance); } } template void checkResults(presult_func_t && presult_func, const Result & results, const DOFs & dofs) { auto presult = presult_func(prescribed_gradient(dofs)); for (auto & result : make_view(results, results.getNbComponent() / dim, dim)) { auto diff = result - presult; auto result_error = diff.template norm() / presult.template norm(); EXPECT_NEAR(0, result_error, result_tolerance); } } template void setLinearDOF(V1 && dof, V2 && coord) { for (UInt i = 0; i < dof.size(); ++i) { dof(i) = this->alpha(i, 0); for (UInt j = 0; j < coord.size(); ++j) { dof(i) += this->alpha(i, j + 1) * coord(j); } } } template void checkDOFs(V && dofs) { const auto & coordinates = mesh->getNodes(); Vector ref_dof(dofs.getNbComponent()); for (auto && tuple : zip(make_view(coordinates, dim), make_view(dofs, dofs.getNbComponent()))) { setLinearDOF(ref_dof, std::get<0>(tuple)); auto diff = std::get<1>(tuple) - ref_dof; auto dofs_error = diff.template norm(); EXPECT_NEAR(0, dofs_error, dofs_tolerance); } } protected: std::unique_ptr mesh; std::unique_ptr model; Matrix alpha{{0.01, 0.02, 0.03, 0.04}, {0.05, 0.06, 0.07, 0.08}, {0.09, 0.10, 0.11, 0.12}}; Real gradient_tolerance{1e-13}; Real result_tolerance{1e-13}; Real dofs_tolerance{1e-15}; }; template constexpr ElementType TestPatchTestLinear::type; template constexpr size_t TestPatchTestLinear::dim; #endif /* __AKANTU_PATCH_TEST_LINEAR_FIXTURE_HH__ */