diff --git a/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.cc b/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.cc index 99985a041..72530f60e 100644 --- a/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.cc +++ b/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.cc @@ -1,73 +1,76 @@ /** * @file boundary_functions.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief * * @section LICENSE * * Copyright (©) 2015-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 "boundary_functions.hh" #include "communicator.hh" +#include "element_group.hh" +#include "node_group.hh" +#include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Real integrateResidual(const std::string & sub_boundary_name, const SolidMechanicsModel & model, UInt dir) { Real int_res = 0.; const Mesh & mesh = model.getMesh(); const Array & residual = model.getInternalForce(); const ElementGroup & boundary = mesh.getElementGroup(sub_boundary_name); - for (auto & node : boundary.getNodes()) { + for (auto & node : boundary.getNodeGroup().getNodes()) { bool is_local_node = mesh.isLocalOrMasterNode(node); if (is_local_node) { int_res += residual(node, dir); } } mesh.getCommunicator().allReduce(int_res, SynchronizerOperation::_sum); return int_res; } /* -------------------------------------------------------------------------- */ void boundaryFix(Mesh & mesh, const std::vector & sub_boundary_names) { std::vector::const_iterator it = sub_boundary_names.begin(); std::vector::const_iterator end = sub_boundary_names.end(); for (; it != end; ++it) { if (mesh.element_group_find(*it) == mesh.element_group_end()) { - mesh.createElementGroup( - *it, mesh.getSpatialDimension() - 1); // empty element group + mesh.createElementGroup(*it, mesh.getSpatialDimension() - + 1); // empty element group } } } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.hh b/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.hh index a16426a6c..f5413bdd9 100644 --- a/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.hh +++ b/extra_packages/traction-at-split-node-contact/src/functions/boundary_functions.hh @@ -1,45 +1,55 @@ /** * @file boundary_functions.hh * * @author David Simon Kammer * * @date creation: Fri Jan 04 2013 * @date last modification: Fri Feb 23 2018 * * @brief functions for boundaries * * @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 . * */ /* -------------------------------------------------------------------------- */ -// akantu #include "aka_common.hh" -#include "solid_mechanics_model.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_BOUNDARY_FUNCTIONS_HH__ +#define __AKANTU_BOUNDARY_FUNCTIONS_HH__ + +namespace akantu { +class SolidMechanicsModel; +} namespace akantu { Real integrateResidual(const std::string & sub_boundary_name, const SolidMechanicsModel & model, UInt dir); /// this is a fix so that all subboundaries exist on all procs void boundaryFix(Mesh & mesh, const std::vector & sub_boundary_names); } // namespace akantu + +#endif /* __AKANTU_BOUNDARY_FUNCTIONS_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc index 4e4232b23..7a93bd5ed 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc @@ -1,554 +1,554 @@ /** * @file ntn_contact.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of ntn_contact * * @section LICENSE * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_contact.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ NTNContact::NTNContact(SolidMechanicsModel & model, const ID & id, const MemoryID & memory_id) : NTNBaseContact(model, id, memory_id), masters(0, 1, 0, id + ":masters", std::numeric_limits::quiet_NaN(), "masters"), lumped_boundary_masters(0, 1, 0, id + ":lumped_boundary_masters", std::numeric_limits::quiet_NaN(), "lumped_boundary_masters"), master_elements("master_elements", id, memory_id) { AKANTU_DEBUG_IN(); const Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = this->model.getSpatialDimension(); this->master_elements.initialize(mesh, _nb_component = 1, _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::pairInterfaceNodes(const ElementGroup & slave_boundary, const ElementGroup & master_boundary, UInt surface_normal_dir, const Mesh & mesh, Array & pairs) { AKANTU_DEBUG_IN(); pairs.resize(0); AKANTU_DEBUG_ASSERT(pairs.getNbComponent() == 2, "Array of node pairs should have nb_component = 2," << " but has nb_component = " << pairs.getNbComponent()); UInt dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(surface_normal_dir < dim, "Mesh is of " << dim << " dimensions" << " and cannot have direction " << surface_normal_dir << " for surface normal"); // offset for projection computation Vector offset(dim - 1); for (UInt i = 0, j = 0; i < dim; ++i) { if (surface_normal_dir != i) { offset(j) = i; ++j; } } // find projected node coordinates const Array & coordinates = mesh.getNodes(); // find slave nodes Array proj_slave_coord(slave_boundary.getNbNodes(), dim - 1, 0.); Array slave_nodes(slave_boundary.getNbNodes()); UInt n(0); - for (auto && slave_node : slave_boundary.getNodes()) { + for (auto && slave_node : slave_boundary.getNodeGroup().getNodes()) { for (UInt d = 0; d < dim - 1; ++d) { proj_slave_coord(n, d) = coordinates(slave_node, offset[d]); slave_nodes(n) = slave_node; } ++n; } // find master nodes Array proj_master_coord(master_boundary.getNbNodes(), dim - 1, 0.); Array master_nodes(master_boundary.getNbNodes()); n = 0; - for (auto && master_node : master_boundary.getNodes()) { + for (auto && master_node : master_boundary.getNodeGroup().getNodes()) { for (UInt d = 0; d < dim - 1; ++d) { proj_master_coord(n, d) = coordinates(master_node, offset[d]); master_nodes(n) = master_node; } ++n; } // find minimum distance between slave nodes to define tolerance Real min_dist = std::numeric_limits::max(); for (UInt i = 0; i < proj_slave_coord.size(); ++i) { for (UInt j = i + 1; j < proj_slave_coord.size(); ++j) { Real dist = 0.; for (UInt d = 0; d < dim - 1; ++d) { dist += (proj_slave_coord(i, d) - proj_slave_coord(j, d)) * (proj_slave_coord(i, d) - proj_slave_coord(j, d)); } if (dist < min_dist) { min_dist = dist; } } } min_dist = std::sqrt(min_dist); Real local_tol = 0.1 * min_dist; // find master slave node pairs for (UInt i = 0; i < proj_slave_coord.size(); ++i) { for (UInt j = 0; j < proj_master_coord.size(); ++j) { Real dist = 0.; for (UInt d = 0; d < dim - 1; ++d) { dist += (proj_slave_coord(i, d) - proj_master_coord(j, d)) * (proj_slave_coord(i, d) - proj_master_coord(j, d)); } dist = std::sqrt(dist); if (dist < local_tol) { // it is a pair Vector pair(2); pair[0] = slave_nodes(i); pair[1] = master_nodes(j); pairs.push_back(pair); continue; // found master do not need to search further for this slave } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addSurfacePair(const ID & slave, const ID & master, UInt surface_normal_dir) { AKANTU_DEBUG_IN(); const Mesh & mesh = this->model.getMesh(); const ElementGroup & slave_boundary = mesh.getElementGroup(slave); const ElementGroup & master_boundary = mesh.getElementGroup(master); this->contact_surfaces.insert(&slave_boundary); this->contact_surfaces.insert(&master_boundary); Array pairs(0, 2); NTNContact::pairInterfaceNodes(slave_boundary, master_boundary, surface_normal_dir, this->model.getMesh(), pairs); // eliminate pairs which contain a pbc slave node Array pairs_no_PBC_slaves(0, 2); Array::const_vector_iterator it = pairs.begin(2); Array::const_vector_iterator end = pairs.end(2); for (; it != end; ++it) { const Vector & pair = *it; if (not mesh.isPeriodicSlave(pair(0)) and not mesh.isPeriodicSlave(pair(1))) { pairs_no_PBC_slaves.push_back(pair); } } this->addNodePairs(pairs_no_PBC_slaves); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addNodePairs(const Array & pairs) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(pairs.getNbComponent() == 2, "Array of node pairs should have nb_component = 2," << " but has nb_component = " << pairs.getNbComponent()); UInt nb_pairs = pairs.size(); for (UInt n = 0; n < nb_pairs; ++n) { this->addSplitNode(pairs(n, 0), pairs(n, 1)); } // synchronize with depending nodes findBoundaryElements(this->slaves.getArray(), this->slave_elements); findBoundaryElements(this->masters.getArray(), this->master_elements); updateInternalData(); syncArrays(_added); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::getNodePairs(Array & pairs) const { AKANTU_DEBUG_IN(); pairs.resize(0); AKANTU_DEBUG_ASSERT(pairs.getNbComponent() == 2, "Array of node pairs should have nb_component = 2," << " but has nb_component = " << pairs.getNbComponent()); UInt nb_pairs = this->getNbContactNodes(); for (UInt n = 0; n < nb_pairs; ++n) { Vector pair{this->slaves(n), this->masters(n)}; pairs.push_back(pair); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addSplitNode(UInt slave, UInt master) { AKANTU_DEBUG_IN(); NTNBaseContact::addSplitNode(slave); this->masters.push_back(master); this->lumped_boundary_masters.push_back( std::numeric_limits::quiet_NaN()); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* This function only works for surface elements with one quad point. For surface elements with more quad points, it computes still, but the result might not be what you are looking for. */ void NTNContact::updateNormals() { AKANTU_DEBUG_IN(); // set normals to zero this->normals.clear(); // contact information UInt dim = this->model.getSpatialDimension(); UInt nb_contact_nodes = this->getNbContactNodes(); this->synch_registry->synchronize(_gst_cf_nodal); // synchronize current pos const Array & cur_pos = this->model.getCurrentPosition(); FEEngine & boundary_fem = this->model.getFEEngineBoundary(); const Mesh & mesh = this->model.getMesh(); for (auto ghost_type: ghost_types) { for (auto & type : mesh.elementTypes(dim - 1, ghost_type)) { // compute the normals Array quad_normals(0, dim); boundary_fem.computeNormalsOnIntegrationPoints(cur_pos, quad_normals, type, ghost_type); UInt nb_quad_points = boundary_fem.getNbIntegrationPoints(type, ghost_type); // new version: compute normals only based on master elements (and not all // boundary elements) // ------------------------------------------------------------------------------------- UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type, ghost_type); // loop over contact nodes for (auto & element : (this->master_elements)(type, ghost_type)) { for (UInt q = 0; q < nb_nodes_per_element; ++q) { UInt node = connectivity(element, q); UInt node_index = this->masters.find(node); AKANTU_DEBUG_ASSERT(node_index != UInt(-1), "Could not find node " << node << " in the array!"); for (UInt q = 0; q < nb_quad_points; ++q) { // add quad normal to master normal for (UInt d = 0; d < dim; ++d) { this->normals(node_index, d) += quad_normals(element * nb_quad_points + q, d); } } } } } } Real * master_normals = this->normals.storage(); for (UInt n = 0; n < nb_contact_nodes; ++n) { if (dim == 2) Math::normalize2(&(master_normals[n * dim])); else if (dim == 3) Math::normalize3(&(master_normals[n * dim])); } // // normalize normals // auto nit = this->normals.begin(); // auto nend = this->normals.end(); // for (; nit != nend; ++nit) { // nit->normalize(); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); NTNBaseContact::dumpRestart(file_name); this->masters.dumpRestartFile(file_name); this->lumped_boundary_masters.dumpRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); NTNBaseContact::readRestart(file_name); this->masters.readRestartFile(file_name); this->lumped_boundary_masters.readRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::updateImpedance() { AKANTU_DEBUG_IN(); UInt nb_contact_nodes = getNbContactNodes(); Real delta_t = this->model.getTimeStep(); AKANTU_DEBUG_ASSERT(delta_t != NAN, "Time step is NAN. Have you set it already?"); const Array & mass = this->model.getMass(); for (UInt n = 0; n < nb_contact_nodes; ++n) { UInt master = this->masters(n); UInt slave = this->slaves(n); Real imp = (this->lumped_boundary_masters(n) / mass(master)) + (this->lumped_boundary_slaves(n) / mass(slave)); imp = 2 / delta_t / imp; this->impedance(n) = imp; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::updateLumpedBoundary() { AKANTU_DEBUG_IN(); internalUpdateLumpedBoundary(this->slaves.getArray(), this->slave_elements, this->lumped_boundary_slaves); internalUpdateLumpedBoundary(this->masters.getArray(), this->master_elements, this->lumped_boundary_masters); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::applyContactPressure() { AKANTU_DEBUG_IN(); UInt nb_ntn_pairs = getNbContactNodes(); UInt dim = this->model.getSpatialDimension(); Array & residual = this->model.getInternalForce(); for (UInt n = 0; n < nb_ntn_pairs; ++n) { UInt master = this->masters(n); UInt slave = this->slaves(n); for (UInt d = 0; d < dim; ++d) { residual(master, d) += this->lumped_boundary_masters(n) * this->contact_pressure(n, d); residual(slave, d) -= this->lumped_boundary_slaves(n) * this->contact_pressure(n, d); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::computeRelativeTangentialField( const Array & field, Array & rel_tang_field) const { AKANTU_DEBUG_IN(); // resize arrays to zero rel_tang_field.resize(0); UInt dim = this->model.getSpatialDimension(); auto it_field = field.begin(dim); auto it_normal = this->normals.getArray().begin(dim); Vector rfv(dim); Vector np_rfv(dim); UInt nb_contact_nodes = this->slaves.size(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // nodes UInt slave = this->slaves(n); UInt master = this->masters(n); // relative field vector (slave - master) rfv = Vector(it_field[slave]); rfv -= Vector(it_field[master]); // normal projection of relative field const Vector normal_v = it_normal[n]; np_rfv = normal_v; np_rfv *= rfv.dot(normal_v); // subract normal projection from relative field to get the tangential // projection rfv -= np_rfv; rel_tang_field.push_back(rfv); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::computeRelativeNormalField( const Array & field, Array & rel_normal_field) const { AKANTU_DEBUG_IN(); // resize arrays to zero rel_normal_field.resize(0); UInt dim = this->model.getSpatialDimension(); // Real * field_p = field.storage(); // Real * normals_p = this->normals.storage(); Array::const_iterator> it_field = field.begin(dim); Array::const_iterator> it_normal = this->normals.getArray().begin(dim); Vector rfv(dim); UInt nb_contact_nodes = this->getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // nodes UInt slave = this->slaves(n); UInt master = this->masters(n); // relative field vector (slave - master) rfv = Vector(it_field[slave]); rfv -= Vector(it_field[master]); // length of normal projection of relative field const Vector normal_v = it_normal[n]; rel_normal_field.push_back(rfv.dot(normal_v)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Int NTNContact::getNodeIndex(UInt node) const { AKANTU_DEBUG_IN(); Int slave_i = NTNBaseContact::getNodeIndex(node); Int master_i = this->masters.find(node); AKANTU_DEBUG_OUT(); return std::max(slave_i, master_i); } /* -------------------------------------------------------------------------- */ void NTNContact::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNContact [" << std::endl; NTNBaseContact::printself(stream, indent); stream << space << " + masters : " << std::endl; this->masters.printself(stream, indent + 2); stream << space << " + lumped_boundary_mastres : " << std::endl; this->lumped_boundary_masters.printself(stream, indent + 2); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::syncArrays(SyncChoice sync_choice) { AKANTU_DEBUG_IN(); NTNBaseContact::syncArrays(sync_choice); this->masters.syncElements(sync_choice); this->lumped_boundary_masters.syncElements(sync_choice); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); /* #ifdef AKANTU_USE_IOHELPER const Array & nodal_filter = this->slaves.getArray(); #define ADD_FIELD(field_id, field, type) \ internalAddDumpFieldToDumper(dumper_name, \ field_id, \ new DumperIOHelper::NodalField< type, true, \ Array, \ Array >(field, 0, 0, &nodal_filter)) */ if (field_id == "lumped_boundary_master") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->lumped_boundary_masters.getArray())); } else { NTNBaseContact::addDumpFieldToDumper(dumper_name, field_id); } /* #undef ADD_FIELD #endif */ AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntrf_contact.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntrf_contact.cc index f0b0c5dea..c63898116 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntrf_contact.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntrf_contact.cc @@ -1,322 +1,322 @@ /** * @file ntrf_contact.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief * * @section LICENSE * * Copyright (©) 2015-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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntrf_contact.hh" namespace akantu { /* -------------------------------------------------------------------------- */ NTRFContact::NTRFContact(SolidMechanicsModel & model, const ID & id, const MemoryID & memory_id) : NTNBaseContact(model, id, memory_id), reference_point(model.getSpatialDimension()), normal(model.getSpatialDimension()) { AKANTU_DEBUG_IN(); is_ntn_contact = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::setReferencePoint(Real x, Real y, Real z) { AKANTU_DEBUG_IN(); Real coord[3]; coord[0] = x; coord[1] = y; coord[2] = z; UInt dim = this->model.getSpatialDimension(); for (UInt d = 0; d < dim; ++d) this->reference_point(d) = coord[d]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::setNormal(Real x, Real y, Real z) { AKANTU_DEBUG_IN(); UInt dim = this->model.getSpatialDimension(); Real coord[3]; coord[0] = x; coord[1] = y; coord[2] = z; for (UInt d = 0; d < dim; ++d) this->normal(d) = coord[d]; this->normal.normalize(); this->updateNormals(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::addSurface(const ID & surf) { AKANTU_DEBUG_IN(); const Mesh & mesh_ref = this->model.getMesh(); try { const ElementGroup & boundary = mesh_ref.getElementGroup(surf); this->contact_surfaces.insert(&boundary); // find slave nodes - for (auto && node : boundary.getNodes()) { + for (auto && node : boundary.getNodeGroup().getNodes()) { if (not mesh_ref.isPeriodicSlave(node)) { this->addSplitNode(node); } } } catch (debug::Exception & e) { AKANTU_DEBUG_INFO("NTRFContact addSurface did not found subboundary " << surf << " and ignored it. Other procs might have it :)"); } // synchronize with depending nodes findBoundaryElements(this->slaves.getArray(), this->slave_elements); updateInternalData(); syncArrays(_added); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::addNodes(Array & nodes) { AKANTU_DEBUG_IN(); UInt nb_nodes = nodes.size(); UInt nb_compo = nodes.getNbComponent(); for (UInt n = 0; n < nb_nodes; ++n) { for (UInt c = 0; c < nb_compo; ++c) { this->addSplitNode(nodes(n, c)); } } // synchronize with depending nodes findBoundaryElements(this->slaves.getArray(), this->slave_elements); updateInternalData(); syncArrays(_added); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::updateNormals() { AKANTU_DEBUG_IN(); // normal is the same for all slaves this->normals.set(this->normal); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::updateImpedance() { AKANTU_DEBUG_IN(); UInt nb_contact_nodes = getNbContactNodes(); Real delta_t = this->model.getTimeStep(); AKANTU_DEBUG_ASSERT(delta_t != NAN, "Time step is NAN. Have you set it already?"); const Array & mass = this->model.getMass(); for (UInt n = 0; n < nb_contact_nodes; ++n) { UInt slave = this->slaves(n); Real imp = this->lumped_boundary_slaves(n) / mass(slave); imp = 2 / delta_t / imp; this->impedance(n) = imp; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::computeRelativeTangentialField( const Array & field, Array & rel_tang_field) const { AKANTU_DEBUG_IN(); // resize arrays to zero rel_tang_field.resize(0); UInt dim = this->model.getSpatialDimension(); Array::const_iterator> it_field = field.begin(dim); Array::const_iterator> it_normal = this->normals.getArray().begin(dim); Vector rfv(dim); Vector np_rfv(dim); UInt nb_contact_nodes = this->slaves.size(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // nodes UInt node = this->slaves(n); // relative field vector rfv = it_field[node]; ; // normal projection of relative field const Vector & normal_v = it_normal[n]; np_rfv = normal_v; np_rfv *= rfv.dot(normal_v); // subtract normal projection from relative field to get the tangential // projection rfv -= np_rfv; rel_tang_field.push_back(rfv); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::computeNormalGap(Array & gap) const { AKANTU_DEBUG_IN(); gap.resize(0); UInt dim = this->model.getSpatialDimension(); Array::const_iterator> it_cur_pos = this->model.getCurrentPosition().begin(dim); Array::const_iterator> it_normal = this->normals.getArray().begin(dim); Vector gap_v(dim); UInt nb_contact_nodes = this->getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // nodes UInt node = this->slaves(n); // gap vector gap_v = it_cur_pos[node]; gap_v -= this->reference_point; // length of normal projection of gap vector const Vector & normal_v = it_normal[n]; gap.push_back(gap_v.dot(normal_v)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::computeRelativeNormalField( const Array & field, Array & rel_normal_field) const { AKANTU_DEBUG_IN(); // resize arrays to zero rel_normal_field.resize(0); UInt dim = this->model.getSpatialDimension(); Array::const_iterator> it_field = field.begin(dim); Array::const_iterator> it_normal = this->normals.getArray().begin(dim); UInt nb_contact_nodes = this->getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // nodes UInt node = this->slaves(n); const Vector & field_v = it_field[node]; const Vector & normal_v = it_normal[n]; rel_normal_field.push_back(field_v.dot(normal_v)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTRFContact [" << std::endl; NTNBaseContact::printself(stream, indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTRFContact::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); /* #ifdef AKANTU_USE_IOHELPER const Array & nodal_filter = this->slaves.getArray(); #define ADD_FIELD(field_id, field, type) \ internalAddDumpFieldToDumper(dumper_name, \ field_id, \ new DumperIOHelper::NodalField< type, true, \ Array, \ Array >(field, 0, 0, &nodal_filter)) */ /* if(field_id == "displacement") { ADD_FIELD(field_id, this->model.getDisplacement(), Real); } else if(field_id == "contact_pressure") { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::NodalField(this->contact_pressure.getArray())); } else {*/ NTNBaseContact::addDumpFieldToDumper(dumper_name, field_id); //} /* #undef ADD_FIELD #endif */ AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/io/mesh_io.cc b/src/io/mesh_io.cc index 0ab01ef38..ab1429c04 100644 --- a/src/io/mesh_io.cc +++ b/src/io/mesh_io.cc @@ -1,141 +1,139 @@ /** * @file mesh_io.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Feb 01 2018 * * @brief common part for all mesh io 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 "mesh_io.hh" #include "aka_common.hh" #include "aka_iterators.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MeshIO::MeshIO() { canReadSurface = false; canReadExtendedData = false; } /* -------------------------------------------------------------------------- */ MeshIO::~MeshIO() = default; /* -------------------------------------------------------------------------- */ std::unique_ptr MeshIO::getMeshIO(const std::string & filename, const MeshIOType & type) { MeshIOType t = type; if (type == _miot_auto) { std::string::size_type idx = filename.rfind('.'); std::string ext; if (idx != std::string::npos) { ext = filename.substr(idx + 1); } if (ext == "msh") { t = _miot_gmsh; } else if (ext == "diana") { t = _miot_diana; } else AKANTU_EXCEPTION("Cannot guess the type of file of " << filename << " (ext " << ext << "). " << "Please provide the MeshIOType to the read function"); } switch (t) { case _miot_gmsh: return std::make_unique(); #if defined(AKANTU_STRUCTURAL_MECHANICS) case _miot_gmsh_struct: return std::make_unique(); #endif case _miot_diana: return std::make_unique(); default: return nullptr; } } /* -------------------------------------------------------------------------- */ void MeshIO::read(const std::string & filename, Mesh & mesh, const MeshIOType & type) { std::unique_ptr mesh_io = getMeshIO(filename, type); mesh_io->read(filename, mesh); } /* -------------------------------------------------------------------------- */ void MeshIO::write(const std::string & filename, Mesh & mesh, const MeshIOType & type) { std::unique_ptr mesh_io = getMeshIO(filename, type); mesh_io->write(filename, mesh); } /* -------------------------------------------------------------------------- */ void MeshIO::constructPhysicalNames(const std::string & tag_name, Mesh & mesh) { - if (!phys_name_map.empty()) { + if (!physical_names.empty()) { for (auto type : mesh.elementTypes()) { auto & name_vec = mesh.getDataPointer("physical_names", type); const auto & tags_vec = mesh.getData(tag_name, type); for (auto pair : zip(tags_vec, name_vec)) { auto tag = std::get<0>(pair); auto & name = std::get<1>(pair); - auto map_it = phys_name_map.find(tag); + auto map_it = physical_names.find(tag); - if (map_it == phys_name_map.end()) { + if (map_it == physical_names.end()) { std::stringstream sstm; sstm << tag; name = sstm.str(); } else { name = map_it->second; } } } } } /* -------------------------------------------------------------------------- */ void MeshIO::printself(std::ostream & stream, int indent) const { - std::string space; - for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) - ; - - if (phys_name_map.size()) { + std::string space(AKANTU_INDENT, indent); + + if (physical_names.size()) { stream << space << "Physical map:" << std::endl; - for (auto & pair : phys_name_map) { + for (auto & pair : physical_names) { stream << space << pair.first << ": " << pair.second << std::endl; } } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/io/mesh_io.hh b/src/io/mesh_io.hh index 03ceb12fe..169e69afd 100644 --- a/src/io/mesh_io.hh +++ b/src/io/mesh_io.hh @@ -1,116 +1,116 @@ /** * @file mesh_io.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Aug 09 2017 * * @brief interface of a mesh io class, reader and writer * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_IO_HH__ #define __AKANTU_MESH_IO_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_accessor.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class MeshIO { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshIO(); virtual ~MeshIO(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void read(const std::string & filename, Mesh & mesh, const MeshIOType & type); void write(const std::string & filename, Mesh & mesh, const MeshIOType & type); /// read a mesh from the file virtual void read(__attribute__((unused)) const std::string & filename, __attribute__((unused)) Mesh & mesh) {} /// write a mesh to a file virtual void write(__attribute__((unused)) const std::string & filename, __attribute__((unused)) const Mesh & mesh) {} /// function to request the manual construction of the physical names maps virtual void constructPhysicalNames(const std::string & tag_name, Mesh & mesh); /// method to permit to be printed to a generic stream virtual void printself(std::ostream & stream, int indent = 0) const; /// static contruction of a meshio object static std::unique_ptr getMeshIO(const std::string & filename, const MeshIOType & type); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - std::map & getPhysicalNameMap() { return phys_name_map; } + auto & getPhysicalNames() { return this->physical_names; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: bool canReadSurface; bool canReadExtendedData; /// correspondance between a tag and physical names (if applicable) - std::map phys_name_map; + std::map physical_names; }; /* -------------------------------------------------------------------------- */ inline std::ostream & operator<<(std::ostream & stream, const MeshIO & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ } // namespace akantu #include "mesh_io_diana.hh" #include "mesh_io_msh.hh" #if defined(AKANTU_STRUCTURAL_MECHANICS) #include "mesh_io_msh_struct.hh" #endif #endif /* __AKANTU_MESH_IO_HH__ */ diff --git a/src/io/mesh_io/mesh_io_msh.cc b/src/io/mesh_io/mesh_io_msh.cc index 40de37cc6..03cd99c56 100644 --- a/src/io/mesh_io/mesh_io_msh.cc +++ b/src/io/mesh_io/mesh_io_msh.cc @@ -1,731 +1,916 @@ /** * @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 - -/* -------------------------------------------------------------------------- */ +#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; /* -------------------------------------------------------------------------- */ -static void my_getline(std::ifstream & infile, std::string & str) { - std::string tmp_str; - std::getline(infile, tmp_str); - str = trim(tmp_str); -} - -/* -------------------------------------------------------------------------- */ -void MeshIOMSH::read(const std::string & filename, Mesh & mesh) { - - MeshAccessor mesh_accessor(mesh); +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; + + Mesh & mesh; + MeshAccessor mesh_accessor; + + std::multimap 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); + } + } - std::ifstream infile; - infile.open(filename.c_str()); + ~File() { infile.close(); } - std::string line; - UInt first_node_number = std::numeric_limits::max(), - last_node_number = 0, file_format = 1, current_line = 0; - bool has_physical_names = false; + auto good() { return infile.good(); } - if (!infile.good()) { - AKANTU_EXCEPTION("Cannot open file " << filename); - } + 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; - std::map> readers; + return std::stringstream(line); + } - readers["$MeshFormat"] = [&](const std::string &) { - my_getline(infile, line); /// the format line - std::stringstream sstr(line); - int version; - sstr >> version; - if (version > 2) - AKANTU_ERROR("This reader can not read version " - << version << " of the MSH file format"); - Int format; - sstr >> format; - if (format != 0) - AKANTU_ERROR("This reader can only read ASCII files."); - my_getline(infile, line); /// the end of block line - current_line += 2; - file_format = 2; + 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); - my_getline(infile, line); - std::stringstream sstr(line); - sstr >> nb_nodes; - current_line++; - - Array & nodes = mesh_accessor.getNodes(); + Array & nodes = file.mesh_accessor.getNodes(); nodes.resize(nb_nodes); - mesh_accessor.setNbGlobalNodes(nb_nodes); + file.mesh_accessor.setNbGlobalNodes(nb_nodes); - UInt index; - Real coord[3]; + size_t index; + double coord[3]; UInt spatial_dimension = nodes.getNbComponent(); /// for each node, read the coordinates for (UInt i = 0; i < nb_nodes; ++i) { UInt offset = i * spatial_dimension; - my_getline(infile, line); - std::stringstream sstr_node(line); - sstr_node >> index >> coord[0] >> coord[1] >> coord[2]; - current_line++; + file.read_line(index, coord[0], coord[1], coord[2]); - first_node_number = std::min(first_node_number, index); - last_node_number = std::max(last_node_number, index); + file.first_node_number = std::min(file.first_node_number, index); + file.last_node_number = std::max(file.last_node_number, index); /// read the coordinates for (UInt j = 0; j < spatial_dimension; ++j) nodes.storage()[offset + j] = coord[j]; } - my_getline(infile, line); /// the end of block line }; readers["$ELM"] = readers["$Elements"] = [&](const std::string &) { UInt nb_elements; - - std::vector read_order; - - my_getline(infile, line); - std::stringstream sstr(line); - sstr >> nb_elements; - current_line++; + file.read_line(nb_elements); Int index; UInt msh_type; - ElementType akantu_type, akantu_type_old = _not_defined; - Array * connectivity = nullptr; - UInt node_per_element = 0; + ElementType akantu_type; for (UInt i = 0; i < nb_elements; ++i) { - my_getline(infile, line); - std::stringstream sstr_elem(line); - current_line++; + 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]; + this->_msh_to_akantu_element_types[MSHElementType(msh_type)]; if (akantu_type == _not_defined) { AKANTU_DEBUG_WARNING("Unsuported element kind " - << msh_type << " at line " << current_line); + << msh_type << " at line " << file.current_line); continue; } - if (akantu_type != akantu_type_old) { - connectivity = &mesh_accessor.getConnectivity(akantu_type); - // connectivity->resize(0); - - node_per_element = connectivity->getNbComponent(); - akantu_type_old = akantu_type; - read_order = this->_read_order[akantu_type]; - } + 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_format == 2) { + if (file.version == 2) { UInt nb_tags; sstr_elem >> nb_tags; for (UInt j = 0; j < nb_tags; ++j) { Int tag; sstr_elem >> tag; - std::stringstream sstr_tag_name; - sstr_tag_name << "tag_" << j; - Array & data = mesh.getDataPointer( - sstr_tag_name.str(), akantu_type, _not_ghost); + + auto & data = file.mesh_accessor.template getData( + "tag_" + std::to_string(j), akantu_type); data.push_back(tag); } - } else if (file_format == 1) { - Int tag; - sstr_elem >> tag; // reg-phys - std::string tag_name = "tag_0"; - Array * data = - &mesh.getDataPointer(tag_name, akantu_type, _not_ghost); - data->push_back(tag); - - sstr_elem >> tag; // reg-elem - tag_name = "tag_1"; - data = &mesh.getDataPointer(tag_name, akantu_type, _not_ghost); - data->push_back(tag); - - sstr_elem >> tag; // number-of-nodes + } else if (file.version == 1) { + 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); } 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 <= last_node_number, - "Node number not in range : line " << current_line); + AKANTU_DEBUG_ASSERT(node_index <= file.last_node_number, + "Node number not in range : line " + << file.current_line); - node_index -= first_node_number; + node_index -= file.first_node_number; local_connect(read_order[j]) = node_index; } - connectivity->push_back(local_connect); + connectivity.push_back(local_connect); } - my_getline(infile, line); /// the end of block line }; - readers["$PhysicalNames"] = [&](const std::string &) { - has_physical_names = true; - my_getline(infile, line); /// the format line - std::stringstream sstr(line); + readers["$NodeData"] = [&](const std::string &) { + /* $NodeData + number-of-string-tags + < "string-tag" > + … + number-of-real-tags + < real-tag > + … + number-of-integer-tags + < integer-tag > + … + node-number value … + … + $EndNodeData */ - UInt num_of_phys_names; - sstr >> num_of_phys_names; + auto read_data_tags = [&](auto x) { + UInt number_of_tags{0}; + file.read_line(number_of_tags); + std::vector tags(number_of_tags); - for (UInt k(0); k < num_of_phys_names; k++) { - my_getline(infile, line); - std::stringstream sstr_phys_name(line); - UInt phys_name_id; - UInt phys_dim; + for (auto && tag : tags) { + file.read_line(tag); + } + return tags; + }; - sstr_phys_name >> phys_dim >> phys_name_id; + auto && string_tags = read_data_tags(std::string{}); + auto && real_tags [[gnu::unused]] = read_data_tags(double{}); + auto && int_tags = read_data_tags(int{}); - std::size_t b = line.find('\"'); - std::size_t e = line.rfind('\"'); - std::string phys_name = line.substr(b + 1, e - b - 1); + auto && data = file.mesh.template registerNodalData( + trim(string_tags[0], '"'), int_tags[1]); + data.resize(file.mesh.getNbNodes(), 0.); + for (auto n [[gnu::unused]] : arange(int_tags[2])) { + auto && sstr = file.get_line(); - phys_name_map[phys_name_id] = phys_name; + int node; + double value; + sstr >> node; + for (auto c : arange(int_tags[1])) { + sstr >> value; + data(node - file.first_node_number, c) = value; + } } - my_getline(infile, line); /// the end of block line }; readers["$Periodic"] = [&](const std::string &) { UInt nb_periodic_entities; - my_getline(infile, line); - - std::stringstream sstr(line); - sstr >> nb_periodic_entities; + file.read_line(nb_periodic_entities); - mesh_accessor.getNodesFlags().resize(mesh.getNbNodes(), - NodeFlag::_normal); + file.mesh_accessor.getNodesFlags().resize(file.mesh.getNbNodes(), + NodeFlag::_normal); for (UInt p = 0; p < nb_periodic_entities; ++p) { // dimension slave-tag master-tag - my_getline(infile, line); - UInt dimension; - { - std::stringstream sstr(line); - sstr >> dimension; - } + file.read_line(dimension); // transformation - my_getline(infile, line); + file.get_line(); // nb nodes - my_getline(infile, line); UInt nb_nodes; - { - std::stringstream sstr(line); - sstr >> nb_nodes; - } + file.read_line(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { // slave master - my_getline(infile, line); + auto && sstr = file.get_line(); - // The info in the mesh seem inconsistent so they are ignored for know. + // The info in the mesh seem inconsistent so they are ignored for now. continue; - if (dimension == mesh.getSpatialDimension() - 1) { + if (dimension == file.mesh.getSpatialDimension() - 1) { UInt slave, master; - std::stringstream sstr(line); + sstr >> slave; sstr >> master; - mesh_accessor.addPeriodicSlave(slave, master); + file.mesh_accessor.addPeriodicSlave(slave, master); } } } // mesh_accessor.markMeshPeriodic(); - my_getline(infile, line); }; +} - readers["$NodeData"] = [&](const std::string &) { - /* $NodeData - number-of-string-tags - < "string-tag" > - … - number-of-real-tags - < real-tag > - … - number-of-integer-tags - < integer-tag > - … - node-number value … - … - $EndNodeData */ +/* -------------------------------------------------------------------------- */ +template +void MeshIOMSH::populateReaders4(File & file, Readers & readers) { + static std::map entity_type{ + {0, "points"}, + {1, "curve"}, + {2, "surface"}, + {3, "volume"}, + }; - auto read_data_tags = [&](auto && tags) { - my_getline(infile, line); /// number of tags - UInt number_of_tags{0}; - std::stringstream sstr(line); - sstr >> number_of_tags; - tags.resize(number_of_tags); + 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 && tag : tags) { - my_getline(infile, line); - std::stringstream sstr(line); - sstr >> tag; + 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)); + } } - return tags; - }; + } + }; - auto && string_tags = read_data_tags(std::vector()); - auto && real_tags[[gnu::unused]] = read_data_tags(std::vector()); - auto && int_tags = read_data_tags(std::vector()); + 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"); + } - auto && data = mesh.registerNodalData(trim(string_tags[0], '"'), int_tags[1]); - data.resize(mesh.getNbNodes(), 0.); - for (auto n [[gnu::unused]] : arange(int_tags[2])) { - my_getline(infile, line); - std::stringstream sstr(line); - int node; - double value; - sstr >> node; - for (auto c : arange(int_tags[1])) { - sstr >> value; - data(node - first_node_number, c) = value; + 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)); + + 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); + } + } 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); + 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]; + + 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); + for (auto it = range.first; it != range.second; ++it) { + file.mesh.getElementGroup(this->physical_names[it->second]).add(elem); + } } } + }; +} + +/* -------------------------------------------------------------------------- */ +void MeshIOMSH::read(const std::string & filename, Mesh & mesh) { + + File file(filename, mesh); + + std::map> readers; - my_getline(infile, line); + 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) { + int size_of_size_t; + sstr >> size_of_size_t; + if (size_of_size_t > int(sizeof(UInt))) { + AKANTU_DEBUG_WARNING("The size of the indexes in akantu are to small " + "to read this file (akantu " + << sizeof(UInt) << " vs. msh file " + << size_of_size_t << ")"); + } + } + + if (file.version < 4) { + this->populateReaders2(file, readers); + } else { + this->populateReaders4(file, readers); + } }; - readers["Unsupported"] = [&](const std::string & block) { - AKANTU_DEBUG_WARNING("Unsuported block_kind " << line << " at line " - << current_line); + 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 (line != end_block) { - my_getline(infile, line); - current_line++; + while (file.line != end_block) { + file.get_line(); } }; - while (infile.good()) { - my_getline(infile, line); - current_line++; + while (file.good()) { + std::string block; + file.read_line(block); - auto && it = readers.find(line); + auto && it = readers.find(block); if (it != readers.end()) { - it->second(line); - } else if(line.size() != 0) { - readers["Unsupported"](line); + 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); - infile.close(); - - this->constructPhysicalNames("tag_0", mesh); - - if (has_physical_names) - mesh.createGroupsFromMeshData("physical_names"); + if (file.version <= 2) { + 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.1 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"; UInt element_idx = 1; for (auto && type : mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { const auto & connectivity = mesh.getConnectivity(type, _not_ghost); 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 (UInt i = 0; i < connectivity.size(); ++i) { UInt offset = i * connectivity.getNbComponent(); 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][i]; else outfile << " 0"; for (UInt j = 0; j < connectivity.getNbComponent(); ++j) { outfile << " " << connectivity.storage()[offset + j] + 1; } outfile << "\n"; element_idx++; } } outfile << "$EndElements" << "\n"; if (mesh.hasData(MeshDataType::_nodal)) { auto && tags = mesh.getTagNames(); for (auto && tag : tags) { if (mesh.getTypeCode(tag, MeshDataType::_nodal) != _tc_real) 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"; } } outfile.close(); } -/* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- + */ } // namespace akantu diff --git a/src/io/mesh_io/mesh_io_msh.hh b/src/io/mesh_io/mesh_io_msh.hh index 568f5cca6..608e951ef 100644 --- a/src/io/mesh_io/mesh_io_msh.hh +++ b/src/io/mesh_io/mesh_io_msh.hh @@ -1,109 +1,116 @@ /** * @file mesh_io_msh.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Nov 08 2017 * * @brief Read/Write for MSH files * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_IO_MSH_HH__ #define __AKANTU_MESH_IO_MSH_HH__ /* -------------------------------------------------------------------------- */ #include "mesh_io.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class MeshIOMSH : public MeshIO { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshIOMSH(); ~MeshIOMSH() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read a mesh from the file void read(const std::string & filename, Mesh & mesh) override; /// write a mesh to a file void write(const std::string & filename, const Mesh & mesh) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// MSH element types enum MSHElementType { _msh_not_defined = 0, _msh_segment_2 = 1, // 2-node line. _msh_triangle_3 = 2, // 3-node triangle. _msh_quadrangle_4 = 3, // 4-node quadrangle. _msh_tetrahedron_4 = 4, // 4-node tetrahedron. _msh_hexahedron_8 = 5, // 8-node hexahedron. _msh_prism_1 = 6, // 6-node prism. _msh_pyramid_1 = 7, // 5-node pyramid. _msh_segment_3 = 8, // 3-node second order line _msh_triangle_6 = 9, // 6-node second order triangle _msh_quadrangle_9 = 10, // 9-node second order quadrangle _msh_tetrahedron_10 = 11, // 10-node second order tetrahedron _msh_hexahedron_27 = 12, // 27-node second order hexahedron _msh_prism_18 = 13, // 18-node second order prism _msh_pyramid_14 = 14, // 14-node second order pyramid _msh_point = 15, // 1-node point. _msh_quadrangle_8 = 16, // 8-node second order quadrangle _msh_hexahedron_20 = 17, // 20-node second order hexahedron _msh_prism_15 = 18 // 15-node second order prism }; #define MAX_NUMBER_OF_NODE_PER_ELEMENT 10 // tetrahedron of second order /// order in witch element as to be read std::map> _read_order; /// number of nodes per msh element std::map _msh_nodes_per_elem; /// correspondence between msh element types and akantu element types std::map _msh_to_akantu_element_types; /// correspondence between akantu element types and msh element types std::map _akantu_to_msh_element_types; + +protected: + template + void populateReaders2(File & file, Readers & readers); + + template + void populateReaders4(File & file, Readers & readers); }; -} // akantu +} // namespace akantu #endif /* __AKANTU_MESH_IO_MSH_HH__ */ diff --git a/src/mesh/element_group.cc b/src/mesh/element_group.cc index c29d194a9..0870259a4 100644 --- a/src/mesh/element_group.cc +++ b/src/mesh/element_group.cc @@ -1,188 +1,194 @@ /** * @file element_group.cc * * @author Dana Christen * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Mon Jan 22 2018 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_csr.hh" #include "dumpable.hh" #include "dumpable_inline_impl.hh" #include "group_manager.hh" #include "group_manager_inline_impl.cc" #include "mesh.hh" #include "mesh_utils.hh" #include #include #include #include "element_group.hh" #if defined(AKANTU_USE_IOHELPER) #include "dumper_iohelper_paraview.hh" #endif namespace akantu { + /* -------------------------------------------------------------------------- */ ElementGroup::ElementGroup(const std::string & group_name, const Mesh & mesh, NodeGroup & node_group, UInt dimension, const std::string & id, const MemoryID & mem_id) : Memory(id, mem_id), mesh(mesh), name(group_name), elements("elements", id, mem_id), node_group(node_group), dimension(dimension) { AKANTU_DEBUG_IN(); #if defined(AKANTU_USE_IOHELPER) this->registerDumper("paraview_" + group_name, group_name, true); - this->addDumpFilteredMesh(mesh, elements, node_group.getNodes(), dimension); + this->addDumpFilteredMesh(mesh, elements, node_group.getNodes(), _all_dimensions); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementGroup::empty() { elements.free(); } /* -------------------------------------------------------------------------- */ void ElementGroup::append(const ElementGroup & other_group) { AKANTU_DEBUG_IN(); node_group.append(other_group.node_group); /// loop on all element types in all dimensions for (auto ghost_type : ghost_types) { for (auto type : other_group.elementTypes(_ghost_type = ghost_type, _element_kind = _ek_not_defined)) { const Array & other_elem_list = other_group.elements(type, ghost_type); UInt nb_other_elem = other_elem_list.size(); Array * elem_list; UInt nb_elem = 0; /// create current type if doesn't exists, otherwise get information if (elements.exists(type, ghost_type)) { elem_list = &elements(type, ghost_type); nb_elem = elem_list->size(); } else { elem_list = &(elements.alloc(0, 1, type, ghost_type)); } /// append new elements to current list elem_list->resize(nb_elem + nb_other_elem); std::copy(other_elem_list.begin(), other_elem_list.end(), elem_list->begin() + nb_elem); /// remove duplicates std::sort(elem_list->begin(), elem_list->end()); Array::iterator<> end = std::unique(elem_list->begin(), elem_list->end()); elem_list->resize(end - elem_list->begin()); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ElementGroup::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "ElementGroup [" << std::endl; stream << space << " + name: " << name << std::endl; stream << space << " + dimension: " << dimension << std::endl; elements.printself(stream, indent + 1); node_group.printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void ElementGroup::optimize() { // increasing the locality of data when iterating on the element of a group for (auto ghost_type : ghost_types) { for (auto type : elements.elementTypes(_ghost_type = ghost_type)) { Array & els = elements(type, ghost_type); std::sort(els.begin(), els.end()); Array::iterator<> end = std::unique(els.begin(), els.end()); els.resize(end - els.begin()); } } node_group.optimize(); } /* -------------------------------------------------------------------------- */ void ElementGroup::fillFromNodeGroup() { CSR node_to_elem; MeshUtils::buildNode2Elements(this->mesh, node_to_elem, this->dimension); std::set seen; Array::const_iterator<> itn = this->node_group.begin(); Array::const_iterator<> endn = this->node_group.end(); for (; itn != endn; ++itn) { CSR::iterator ite = node_to_elem.begin(*itn); CSR::iterator ende = node_to_elem.end(*itn); for (; ite != ende; ++ite) { const Element & elem = *ite; if (this->dimension != _all_dimensions && this->dimension != Mesh::getSpatialDimension(elem.type)) continue; if (seen.find(elem) != seen.end()) continue; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(elem.type); Array::const_iterator> conn_it = this->mesh.getConnectivity(elem.type, elem.ghost_type) .begin(nb_nodes_per_element); const Vector & conn = conn_it[elem.element]; UInt count = 0; for (UInt n = 0; n < conn.size(); ++n) { count += (this->node_group.getNodes().find(conn(n)) != UInt(-1) ? 1 : 0); } if (count == nb_nodes_per_element) this->add(elem); seen.insert(elem); } } this->optimize(); } +/* -------------------------------------------------------------------------- */ +void ElementGroup::addDimension(UInt dimension) { + this->dimension = std::max(dimension, this->dimension); +} + /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/mesh/element_group.hh b/src/mesh/element_group.hh index 168bc7a12..84a1f2a2b 100644 --- a/src/mesh/element_group.hh +++ b/src/mesh/element_group.hh @@ -1,193 +1,200 @@ /** * @file element_group.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Wed Nov 08 2017 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "dumpable.hh" #include "element_type_map.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_GROUP_HH__ #define __AKANTU_ELEMENT_GROUP_HH__ namespace akantu { class Mesh; class Element; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class ElementGroup : private Memory, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ElementGroup(const std::string & name, const Mesh & mesh, NodeGroup & node_group, UInt dimension = _all_dimensions, const std::string & id = "element_group", const MemoryID & memory_id = 0); /* ------------------------------------------------------------------------ */ /* Type definitions */ /* ------------------------------------------------------------------------ */ public: using ElementList = ElementTypeMapArray; using NodeList = Array; /* ------------------------------------------------------------------------ */ /* Element iterator */ /* ------------------------------------------------------------------------ */ #ifndef SWIG using type_iterator = ElementList::type_iterator; [[deprecated("Use elementTypes instead")]] inline type_iterator firstType(UInt dim = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_regular) const; [[deprecated("Use elementTypes instead")]] inline type_iterator lastType(UInt dim = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_regular) const; template inline decltype(auto) elementTypes(pack &&... _pack) const { return elements.elementTypes(_pack...); } #endif using const_element_iterator = Array::const_iterator; inline const_element_iterator begin(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; inline const_element_iterator end(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// empty the element group void empty(); /// append another group to this group /// BE CAREFUL: it doesn't conserve the element order void append(const ElementGroup & other_group); /// add an element to the group. By default the it does not add the nodes to /// the group inline void add(const Element & el, bool add_nodes = false, bool check_for_duplicate = true); /// \todo fix the default for add_nodes : make it coherent with the other /// method inline void add(const ElementType & type, UInt element, const GhostType & ghost_type = _not_ghost, bool add_nodes = true, bool check_for_duplicate = true); inline void addNode(UInt node_id, bool check_for_duplicate = true); inline void removeNode(UInt node_id); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// fill the elements based on the underlying node group. virtual void fillFromNodeGroup(); // sort and remove duplicated values void optimize(); + /// change the dimension if needed + void addDimension(UInt dimension); + private: inline void addElement(const ElementType & elem_type, UInt elem_id, const GhostType & ghost_type); friend class GroupManager; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: const Array & getElements(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; AKANTU_GET_MACRO(Elements, elements, const ElementTypeMapArray &); - AKANTU_GET_MACRO(Nodes, node_group.getNodes(), const Array &); + AKANTU_GET_MACRO_NOT_CONST(Elements, elements, ElementTypeMapArray &); + +// AKANTU_GET_MACRO(Nodes, node_group.getNodes(), const Array &); + AKANTU_GET_MACRO(NodeGroup, node_group, const NodeGroup &); AKANTU_GET_MACRO_NOT_CONST(NodeGroup, node_group, NodeGroup &); + AKANTU_GET_MACRO(Dimension, dimension, UInt); AKANTU_GET_MACRO(Name, name, std::string); inline UInt getNbNodes() const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Mesh to which this group belongs const Mesh & mesh; /// name of the group std::string name; /// list of elements composing the group ElementList elements; /// sub list of nodes which are composing the elements NodeGroup & node_group; /// group dimension UInt dimension; /// empty arry for the iterator to work when an element type not present Array empty_elements; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const ElementGroup & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "element.hh" #include "element_group_inline_impl.cc" #endif /* __AKANTU_ELEMENT_GROUP_HH__ */ diff --git a/src/mesh/group_manager.cc b/src/mesh/group_manager.cc index 2955d59dd..5c918c92f 100644 --- a/src/mesh/group_manager.cc +++ b/src/mesh/group_manager.cc @@ -1,1005 +1,1063 @@ /** * @file group_manager.cc * * @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: Tue Feb 20 2018 * * @brief Stores information about ElementGroup and NodeGroup * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "group_manager.hh" #include "aka_csr.hh" #include "data_accessor.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "mesh_utils.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ GroupManager::GroupManager(const Mesh & mesh, const ID & id, const MemoryID & mem_id) : id(id), memory_id(mem_id), mesh(mesh) { AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ GroupManager::~GroupManager() { auto eit = element_groups.begin(); auto eend = element_groups.end(); for (; eit != eend; ++eit) delete (eit->second); auto nit = node_groups.begin(); auto nend = node_groups.end(); for (; nit != nend; ++nit) delete (nit->second); } /* -------------------------------------------------------------------------- */ NodeGroup & GroupManager::createNodeGroup(const std::string & group_name, bool replace_group) { AKANTU_DEBUG_IN(); auto it = node_groups.find(group_name); if (it != node_groups.end()) { if (replace_group) { it->second->empty(); AKANTU_DEBUG_OUT(); return *(it->second); } else AKANTU_EXCEPTION( "Trying to create a node group that already exists:" << group_name); } std::stringstream sstr; sstr << this->id << ":" << group_name << "_node_group"; NodeGroup * node_group = new NodeGroup(group_name, mesh, sstr.str(), memory_id); node_groups[group_name] = node_group; AKANTU_DEBUG_OUT(); return *node_group; } /* -------------------------------------------------------------------------- */ template NodeGroup & GroupManager::createFilteredNodeGroup(const std::string & group_name, const NodeGroup & source_node_group, T & filter) { AKANTU_DEBUG_IN(); NodeGroup & node_group = this->createNodeGroup(group_name); node_group.append(source_node_group); if (T::type == FilterFunctor::_node_filter_functor) { node_group.applyNodeFilter(filter); } else { AKANTU_ERROR("ElementFilter cannot be applied to NodeGroup yet." << " Needs to be implemented."); } AKANTU_DEBUG_OUT(); return node_group; } /* -------------------------------------------------------------------------- */ void GroupManager::destroyNodeGroup(const std::string & group_name) { AKANTU_DEBUG_IN(); auto nit = node_groups.find(group_name); auto nend = node_groups.end(); if (nit != nend) { delete (nit->second); node_groups.erase(nit); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::createElementGroup(const std::string & group_name, UInt dimension, bool replace_group) { AKANTU_DEBUG_IN(); NodeGroup & new_node_group = createNodeGroup(group_name + "_nodes", replace_group); auto it = element_groups.find(group_name); if (it != element_groups.end()) { if (replace_group) { - it->second->empty(); - AKANTU_DEBUG_OUT(); - return *(it->second); + destroyElementGroup(group_name, true); } else AKANTU_EXCEPTION("Trying to create a element group that already exists:" << group_name); } - std::stringstream sstr; - sstr << this->id << ":" << group_name << "_element_group"; - ElementGroup * element_group = new ElementGroup( - group_name, mesh, new_node_group, dimension, sstr.str(), memory_id); + group_name, mesh, new_node_group, dimension, + this->id + ":" + group_name + "_element_group", memory_id); - std::stringstream sstr_nodes; - sstr_nodes << group_name << "_nodes"; - - node_groups[sstr_nodes.str()] = &new_node_group; element_groups[group_name] = element_group; AKANTU_DEBUG_OUT(); return *element_group; } /* -------------------------------------------------------------------------- */ void GroupManager::destroyElementGroup(const std::string & group_name, bool destroy_node_group) { AKANTU_DEBUG_IN(); auto eit = element_groups.find(group_name); auto eend = element_groups.end(); if (eit != eend) { if (destroy_node_group) destroyNodeGroup(eit->second->getNodeGroup().getName()); delete (eit->second); element_groups.erase(eit); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::destroyAllElementGroups(bool destroy_node_groups) { AKANTU_DEBUG_IN(); auto eit = element_groups.begin(); auto eend = element_groups.end(); for (; eit != eend; ++eit) { if (destroy_node_groups) destroyNodeGroup(eit->second->getNodeGroup().getName()); delete (eit->second); } element_groups.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group) { AKANTU_DEBUG_IN(); if (element_groups.find(group_name) != element_groups.end()) AKANTU_EXCEPTION( "Trying to create a element group that already exists:" << group_name); ElementGroup * element_group = new ElementGroup(group_name, mesh, node_group, dimension, id + ":" + group_name + "_element_group", memory_id); element_groups[group_name] = element_group; AKANTU_DEBUG_OUT(); return *element_group; } /* -------------------------------------------------------------------------- */ template ElementGroup & GroupManager::createFilteredElementGroup( const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter) { AKANTU_DEBUG_IN(); ElementGroup * element_group = nullptr; if (T::type == FilterFunctor::_node_filter_functor) { NodeGroup & filtered_node_group = this->createFilteredNodeGroup( group_name + "_nodes", node_group, filter); element_group = &(this->createElementGroup(group_name, dimension, filtered_node_group)); } else if (T::type == FilterFunctor::_element_filter_functor) { AKANTU_ERROR( "Cannot handle an ElementFilter yet. Needs to be implemented."); } AKANTU_DEBUG_OUT(); return *element_group; } /* -------------------------------------------------------------------------- */ class ClusterSynchronizer : public DataAccessor { using DistantIDs = std::set>; public: ClusterSynchronizer(GroupManager & group_manager, UInt element_dimension, std::string cluster_name_prefix, ElementTypeMapArray & element_to_fragment, const ElementSynchronizer & element_synchronizer, UInt nb_cluster) : group_manager(group_manager), element_dimension(element_dimension), cluster_name_prefix(std::move(cluster_name_prefix)), element_to_fragment(element_to_fragment), element_synchronizer(element_synchronizer), nb_cluster(nb_cluster) {} UInt synchronize() { Communicator & comm = Communicator::getStaticCommunicator(); UInt rank = comm.whoAmI(); UInt nb_proc = comm.getNbProc(); /// find starting index to renumber local clusters Array nb_cluster_per_proc(nb_proc); nb_cluster_per_proc(rank) = nb_cluster; comm.allGather(nb_cluster_per_proc); starting_index = std::accumulate(nb_cluster_per_proc.begin(), nb_cluster_per_proc.begin() + rank, 0); UInt global_nb_fragment = std::accumulate(nb_cluster_per_proc.begin() + rank, nb_cluster_per_proc.end(), starting_index); /// create the local to distant cluster pairs with neighbors element_synchronizer.synchronizeOnce(*this, _gst_gm_clusters); /// count total number of pairs Array nb_pairs(nb_proc); // This is potentially a bug for more than // 2**31 pairs, but due to a all gatherv after // it must be int to match MPI interfaces nb_pairs(rank) = distant_ids.size(); comm.allGather(nb_pairs); UInt total_nb_pairs = std::accumulate(nb_pairs.begin(), nb_pairs.end(), 0); /// generate pairs global array UInt local_pair_index = std::accumulate(nb_pairs.storage(), nb_pairs.storage() + rank, 0); Array total_pairs(total_nb_pairs, 2); for (auto & ids : distant_ids) { total_pairs(local_pair_index, 0) = ids.first; total_pairs(local_pair_index, 1) = ids.second; ++local_pair_index; } /// communicate pairs to all processors nb_pairs *= 2; comm.allGatherV(total_pairs, nb_pairs); /// renumber clusters /// generate fragment list std::vector> global_clusters; UInt total_nb_cluster = 0; Array is_fragment_in_cluster(global_nb_fragment, 1, false); std::queue fragment_check_list; while (total_pairs.size() != 0) { /// create a new cluster ++total_nb_cluster; global_clusters.resize(total_nb_cluster); std::set & current_cluster = global_clusters[total_nb_cluster - 1]; UInt first_fragment = total_pairs(0, 0); UInt second_fragment = total_pairs(0, 1); total_pairs.erase(0); fragment_check_list.push(first_fragment); fragment_check_list.push(second_fragment); while (!fragment_check_list.empty()) { UInt current_fragment = fragment_check_list.front(); UInt * total_pairs_end = total_pairs.storage() + total_pairs.size() * 2; UInt * fragment_found = std::find(total_pairs.storage(), total_pairs_end, current_fragment); if (fragment_found != total_pairs_end) { UInt position = fragment_found - total_pairs.storage(); UInt pair = position / 2; UInt other_index = (position + 1) % 2; fragment_check_list.push(total_pairs(pair, other_index)); total_pairs.erase(pair); } else { fragment_check_list.pop(); current_cluster.insert(current_fragment); is_fragment_in_cluster(current_fragment) = true; } } } /// add to FragmentToCluster all local fragments for (UInt c = 0; c < global_nb_fragment; ++c) { if (!is_fragment_in_cluster(c)) { ++total_nb_cluster; global_clusters.resize(total_nb_cluster); std::set & current_cluster = global_clusters[total_nb_cluster - 1]; current_cluster.insert(c); } } /// reorganize element groups to match global clusters for (UInt c = 0; c < global_clusters.size(); ++c) { /// create new element group corresponding to current cluster std::stringstream sstr; sstr << cluster_name_prefix << "_" << c; ElementGroup & cluster = group_manager.createElementGroup(sstr.str(), element_dimension, true); auto it = global_clusters[c].begin(); auto end = global_clusters[c].end(); /// append to current element group all fragments that belong to /// the same cluster if they exist for (; it != end; ++it) { Int local_index = *it - starting_index; if (local_index < 0 || local_index >= Int(nb_cluster)) continue; std::stringstream tmp_sstr; tmp_sstr << "tmp_" << cluster_name_prefix << "_" << local_index; auto eg_it = group_manager.element_group_find(tmp_sstr.str()); AKANTU_DEBUG_ASSERT(eg_it != group_manager.element_group_end(), "Temporary fragment \"" << tmp_sstr.str() << "\" not found"); cluster.append(*(eg_it->second)); group_manager.destroyElementGroup(tmp_sstr.str(), true); } } return total_nb_cluster; } private: /// functions for parallel communications inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override { if (tag == _gst_gm_clusters) return elements.size() * sizeof(UInt); return 0; } inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override { if (tag != _gst_gm_clusters) return; Array::const_iterator<> el_it = elements.begin(); Array::const_iterator<> el_end = elements.end(); for (; el_it != el_end; ++el_it) { const Element & el = *el_it; /// for each element pack its global cluster index buffer << element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index; } } inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override { if (tag != _gst_gm_clusters) return; Array::const_iterator<> el_it = elements.begin(); Array::const_iterator<> el_end = elements.end(); for (; el_it != el_end; ++el_it) { UInt distant_cluster; buffer >> distant_cluster; const Element & el = *el_it; UInt local_cluster = element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index; distant_ids.insert(std::make_pair(local_cluster, distant_cluster)); } } private: GroupManager & group_manager; UInt element_dimension; std::string cluster_name_prefix; ElementTypeMapArray & element_to_fragment; const ElementSynchronizer & element_synchronizer; UInt nb_cluster; DistantIDs distant_ids; UInt starting_index; }; /* -------------------------------------------------------------------------- */ /// \todo this function doesn't work in 1D UInt GroupManager::createBoundaryGroupFromGeometry() { UInt spatial_dimension = mesh.getSpatialDimension(); return createClusters(spatial_dimension - 1, "boundary"); } /* -------------------------------------------------------------------------- */ UInt GroupManager::createClusters( UInt element_dimension, Mesh & mesh_facets, std::string cluster_name_prefix, const GroupManager::ClusteringFilter & filter) { return createClusters(element_dimension, std::move(cluster_name_prefix), filter, mesh_facets); } /* -------------------------------------------------------------------------- */ UInt GroupManager::createClusters( UInt element_dimension, std::string cluster_name_prefix, const GroupManager::ClusteringFilter & filter) { std::unique_ptr mesh_facets; if (!mesh_facets && element_dimension > 0) { MeshAccessor mesh_accessor(const_cast(mesh)); mesh_facets = std::make_unique(mesh.getSpatialDimension(), mesh_accessor.getNodesSharedPtr(), "mesh_facets_for_clusters"); mesh_facets->defineMeshParent(mesh); MeshUtils::buildAllFacets(mesh, *mesh_facets, element_dimension, element_dimension - 1); } return createClusters(element_dimension, std::move(cluster_name_prefix), filter, *mesh_facets); } /* -------------------------------------------------------------------------- */ //// \todo if needed element list construction can be optimized by //// templating the filter class UInt GroupManager::createClusters(UInt element_dimension, const std::string & cluster_name_prefix, const GroupManager::ClusteringFilter & filter, Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt nb_proc = mesh.getCommunicator().getNbProc(); std::string tmp_cluster_name_prefix = cluster_name_prefix; ElementTypeMapArray * element_to_fragment = nullptr; if (nb_proc > 1 && mesh.isDistributed()) { element_to_fragment = new ElementTypeMapArray("element_to_fragment", id, memory_id); element_to_fragment->initialize( mesh, _nb_component = 1, _spatial_dimension = element_dimension, _element_kind = _ek_not_defined, _with_nb_element = true); // mesh.initElementTypeMapArray(*element_to_fragment, 1, element_dimension, // false, _ek_not_defined, true); tmp_cluster_name_prefix = "tmp_" + tmp_cluster_name_prefix; } ElementTypeMapArray seen_elements("seen_elements", id, memory_id); seen_elements.initialize(mesh, _spatial_dimension = element_dimension, _element_kind = _ek_not_defined, _with_nb_element = true); // mesh.initElementTypeMapArray(seen_elements, 1, element_dimension, false, // _ek_not_defined, true); for (auto ghost_type : ghost_types) { Element el; el.ghost_type = ghost_type; for (auto type : mesh.elementTypes(_spatial_dimension = element_dimension, _ghost_type = ghost_type, _element_kind = _ek_not_defined)) { el.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array & seen_elements_array = seen_elements(type, ghost_type); for (UInt e = 0; e < nb_element; ++e) { el.element = e; if (!filter(el)) seen_elements_array(e) = true; } } } Array checked_node(mesh.getNbNodes(), 1, false); UInt nb_cluster = 0; for (auto ghost_type : ghost_types) { Element uns_el; uns_el.ghost_type = ghost_type; for (auto type : mesh.elementTypes(_spatial_dimension = element_dimension, _ghost_type = ghost_type, _element_kind = _ek_not_defined)) { uns_el.type = type; Array & seen_elements_vec = seen_elements(uns_el.type, uns_el.ghost_type); for (UInt e = 0; e < seen_elements_vec.size(); ++e) { // skip elements that have been already seen if (seen_elements_vec(e) == true) continue; // set current element uns_el.element = e; seen_elements_vec(e) = true; /// create a new cluster std::stringstream sstr; sstr << tmp_cluster_name_prefix << "_" << nb_cluster; ElementGroup & cluster = createElementGroup(sstr.str(), element_dimension, true); ++nb_cluster; // point element are cluster by themself if (element_dimension == 0) { cluster.add(uns_el); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(uns_el.type); Vector connect = mesh.getConnectivity(uns_el.type, uns_el.ghost_type) .begin(nb_nodes_per_element)[uns_el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { /// add element's nodes to the cluster UInt node = connect[n]; if (!checked_node(node)) { cluster.addNode(node); checked_node(node) = true; } } continue; } std::queue element_to_add; element_to_add.push(uns_el); /// keep looping until current cluster is complete (no more /// connected elements) while (!element_to_add.empty()) { /// take first element and erase it in the queue Element el = element_to_add.front(); element_to_add.pop(); /// if parallel, store cluster index per element if (nb_proc > 1 && mesh.isDistributed()) (*element_to_fragment)(el.type, el.ghost_type)(el.element) = nb_cluster - 1; /// add current element to the cluster cluster.add(el); const Array & element_to_facet = mesh_facets.getSubelementToElement(el.type, el.ghost_type); UInt nb_facet_per_element = element_to_facet.getNbComponent(); for (UInt f = 0; f < nb_facet_per_element; ++f) { const Element & facet = element_to_facet(el.element, f); if (facet == ElementNull) continue; const std::vector & connected_elements = mesh_facets.getElementToSubelement( facet.type, facet.ghost_type)(facet.element); for (UInt elem = 0; elem < connected_elements.size(); ++elem) { const Element & check_el = connected_elements[elem]; // check if this element has to be skipped if (check_el == ElementNull || check_el == el) continue; Array & seen_elements_vec_current = seen_elements(check_el.type, check_el.ghost_type); if (seen_elements_vec_current(check_el.element) == false) { seen_elements_vec_current(check_el.element) = true; element_to_add.push(check_el); } } } UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); Vector connect = mesh.getConnectivity(el.type, el.ghost_type) .begin(nb_nodes_per_element)[el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { /// add element's nodes to the cluster UInt node = connect[n]; if (!checked_node(node)) { cluster.addNode(node, false); checked_node(node) = true; } } } } } } if (nb_proc > 1 && mesh.isDistributed()) { ClusterSynchronizer cluster_synchronizer( *this, element_dimension, cluster_name_prefix, *element_to_fragment, this->mesh.getElementSynchronizer(), nb_cluster); nb_cluster = cluster_synchronizer.synchronize(); delete element_to_fragment; } if (mesh.isDistributed()) this->synchronizeGroupNames(); AKANTU_DEBUG_OUT(); return nb_cluster; } /* -------------------------------------------------------------------------- */ template void GroupManager::createGroupsFromMeshData(const std::string & dataset_name) { std::set group_names; const auto & datas = mesh.getData(dataset_name); std::map group_dim; for (auto ghost_type : ghost_types) { for (auto type : datas.elementTypes(_ghost_type = ghost_type)) { const Array & dataset = datas(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); AKANTU_DEBUG_ASSERT(dataset.size() == nb_element, "Not the same number of elements (" << type << ":" << ghost_type << ") in the map from MeshData (" << dataset.size() << ") " << dataset_name << " and in the mesh (" << nb_element << ")!"); for (UInt e(0); e < nb_element; ++e) { std::stringstream sstr; sstr << dataset(e); std::string gname = sstr.str(); group_names.insert(gname); auto it = group_dim.find(gname); if (it == group_dim.end()) { group_dim[gname] = mesh.getSpatialDimension(type); } else { it->second = std::max(it->second, mesh.getSpatialDimension(type)); } } } } auto git = group_names.begin(); auto gend = group_names.end(); for (; git != gend; ++git) createElementGroup(*git, group_dim[*git]); if (mesh.isDistributed()) this->synchronizeGroupNames(); Element el; for (auto ghost_type : ghost_types) { el.ghost_type = ghost_type; for (auto type : datas.elementTypes(_ghost_type = ghost_type)) { el.type = type; const Array & dataset = datas(type, ghost_type); UInt nb_element = mesh.getNbElement(type, ghost_type); AKANTU_DEBUG_ASSERT(dataset.size() == nb_element, "Not the same number of elements in the map from " "MeshData and in the mesh!"); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(el.type); Array::const_iterator> cit = mesh.getConnectivity(type, ghost_type).begin(nb_nodes_per_element); for (UInt e(0); e < nb_element; ++e, ++cit) { el.element = e; std::stringstream sstr; sstr << dataset(e); ElementGroup & group = getElementGroup(sstr.str()); group.add(el, false, false); const Vector & connect = *cit; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; group.addNode(node, false); } } } } git = group_names.begin(); for (; git != gend; ++git) { getElementGroup(*git).optimize(); } } template void GroupManager::createGroupsFromMeshData( const std::string & dataset_name); template void GroupManager::createGroupsFromMeshData(const std::string & dataset_name); /* -------------------------------------------------------------------------- */ void GroupManager::createElementGroupFromNodeGroup( const std::string & name, const std::string & node_group_name, UInt dimension) { NodeGroup & node_group = getNodeGroup(node_group_name); ElementGroup & group = createElementGroup(name, dimension, node_group); group.fillFromNodeGroup(); } /* -------------------------------------------------------------------------- */ void GroupManager::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "GroupManager [" << std::endl; std::set node_group_seen; for (auto it(element_group_begin()); it != element_group_end(); ++it) { it->second->printself(stream, indent + 1); node_group_seen.insert(it->second->getNodeGroup().getName()); } for (auto it(node_group_begin()); it != node_group_end(); ++it) { if (node_group_seen.find(it->second->getName()) == node_group_seen.end()) it->second->printself(stream, indent + 1); } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ UInt GroupManager::getNbElementGroups(UInt dimension) const { if (dimension == _all_dimensions) return element_groups.size(); auto it = element_groups.begin(); auto end = element_groups.end(); UInt count = 0; for (; it != end; ++it) count += (it->second->getDimension() == dimension); return count; } /* -------------------------------------------------------------------------- */ void GroupManager::checkAndAddGroups(DynamicCommunicationBuffer & buffer) { AKANTU_DEBUG_IN(); UInt nb_node_group; buffer >> nb_node_group; AKANTU_DEBUG_INFO("Received " << nb_node_group << " node group names"); for (UInt ng = 0; ng < nb_node_group; ++ng) { std::string node_group_name; buffer >> node_group_name; if (node_groups.find(node_group_name) == node_groups.end()) { this->createNodeGroup(node_group_name); } AKANTU_DEBUG_INFO("Received node goup name: " << node_group_name); } UInt nb_element_group; buffer >> nb_element_group; AKANTU_DEBUG_INFO("Received " << nb_element_group << " element group names"); for (UInt eg = 0; eg < nb_element_group; ++eg) { std::string element_group_name; buffer >> element_group_name; std::string node_group_name; buffer >> node_group_name; UInt dim; buffer >> dim; AKANTU_DEBUG_INFO("Received element group name: " << element_group_name << " corresponding to a " << Int(dim) << "D group with node group " << node_group_name); NodeGroup & node_group = *node_groups[node_group_name]; if (element_groups.find(element_group_name) == element_groups.end()) { this->createElementGroup(element_group_name, dim, node_group); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::fillBufferWithGroupNames( DynamicCommunicationBuffer & comm_buffer) const { AKANTU_DEBUG_IN(); // packing node group names; UInt nb_groups = this->node_groups.size(); comm_buffer << nb_groups; AKANTU_DEBUG_INFO("Sending " << nb_groups << " node group names"); auto nnames_it = node_groups.begin(); auto nnames_end = node_groups.end(); for (; nnames_it != nnames_end; ++nnames_it) { std::string node_group_name = nnames_it->first; comm_buffer << node_group_name; AKANTU_DEBUG_INFO("Sending node goupe name: " << node_group_name); } // packing element group names with there associated node group name nb_groups = this->element_groups.size(); comm_buffer << nb_groups; AKANTU_DEBUG_INFO("Sending " << nb_groups << " element group names"); auto gnames_it = this->element_groups.begin(); auto gnames_end = this->element_groups.end(); for (; gnames_it != gnames_end; ++gnames_it) { ElementGroup & element_group = *(gnames_it->second); std::string element_group_name = gnames_it->first; std::string node_group_name = element_group.getNodeGroup().getName(); UInt dim = element_group.getDimension(); comm_buffer << element_group_name; comm_buffer << node_group_name; comm_buffer << dim; AKANTU_DEBUG_INFO("Sending element group name: " << element_group_name << " corresponding to a " << Int(dim) << "D group with the node group " << node_group_name); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GroupManager::synchronizeGroupNames() { AKANTU_DEBUG_IN(); const Communicator & comm = mesh.getCommunicator(); Int nb_proc = comm.getNbProc(); Int my_rank = comm.whoAmI(); if (nb_proc == 1) return; if (my_rank == 0) { for (Int p = 1; p < nb_proc; ++p) { DynamicCommunicationBuffer recv_buffer; comm.receive(recv_buffer, p, p); AKANTU_DEBUG_INFO("Got " << printMemorySize(recv_buffer.size()) << " from proc " << p); this->checkAndAddGroups(recv_buffer); } DynamicCommunicationBuffer comm_buffer; this->fillBufferWithGroupNames(comm_buffer); AKANTU_DEBUG_INFO("Initiating broadcast with " << printMemorySize(comm_buffer.size())); comm.broadcast(comm_buffer); } else { DynamicCommunicationBuffer comm_buffer; this->fillBufferWithGroupNames(comm_buffer); AKANTU_DEBUG_INFO("Sending " << printMemorySize(comm_buffer.size()) << " to proc " << 0); comm.send(comm_buffer, 0, my_rank); DynamicCommunicationBuffer recv_buffer; comm.broadcast(recv_buffer); AKANTU_DEBUG_INFO("Receiving broadcast with " << printMemorySize(recv_buffer.size())); this->checkAndAddGroups(recv_buffer); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const ElementGroup & GroupManager::getElementGroup(const std::string & name) const { auto it = element_group_find(name); if (it == element_group_end()) { AKANTU_EXCEPTION("There are no element groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ ElementGroup & GroupManager::getElementGroup(const std::string & name) { auto it = element_group_find(name); if (it == element_group_end()) { AKANTU_EXCEPTION("There are no element groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const NodeGroup & GroupManager::getNodeGroup(const std::string & name) const { auto it = node_group_find(name); if (it == node_group_end()) { AKANTU_EXCEPTION("There are no node groups named " << name << " associated to the group manager: " << id); } return *(it->second); } /* -------------------------------------------------------------------------- */ NodeGroup & GroupManager::getNodeGroup(const std::string & name) { auto it = node_group_find(name); if (it == node_group_end()) { AKANTU_EXCEPTION("There are no node groups named " << name << " associated to the group manager: " << id); } return *(it->second); } +/* -------------------------------------------------------------------------- */ +template +void GroupManager::renameGroup(GroupsType & groups, const std::string & name, + const std::string & new_name) { + auto it = groups.find(name); + if (it == groups.end()) { + AKANTU_EXCEPTION("There are no group named " + << name << " associated to the group manager: " << id); + } + + auto && group_ptr = it->second; + + group_ptr->name = new_name; + + groups.erase(it); + groups[new_name] = group_ptr; +} + +/* -------------------------------------------------------------------------- */ +void GroupManager::renameElementGroup(const std::string & name, + const std::string & new_name) { + renameGroup(element_groups, name, new_name); +} + +/* -------------------------------------------------------------------------- */ +void GroupManager::renameNodeGroup(const std::string & name, + const std::string & new_name) { + renameGroup(node_groups, name, new_name); +} + +/* -------------------------------------------------------------------------- */ +// template +// void GroupManager::renameGroup(GroupsType & groups, const std::string & name, +// const std::string & new_name) { +// auto it = groups.find(name); +// if (it == groups.end()) { +// AKANTU_EXCEPTION("There are no group named " +// << name << " associated to the group manager: " << id); +// } + +// auto && group_ptr = it->second; + +// group_ptr->name = new_name; + +// groups.erase(it); +// groups[new_name] = group_ptr; +// } + +/* -------------------------------------------------------------------------- */ +void GroupManager::copyElementGroup(const std::string & name, + const std::string & new_name) { + const auto & grp = getElementGroup(name); + auto & new_grp = createElementGroup(new_name, grp.getDimension()); + + new_grp.getElements().copy(grp.getElements()); +} + +/* -------------------------------------------------------------------------- */ +void GroupManager::copyNodeGroup(const std::string & name, + const std::string & new_name) { + const auto & grp = getNodeGroup(name); + auto & new_grp = createNodeGroup(new_name); + + new_grp.getNodes().copy(grp.getNodes()); +} + } // namespace akantu diff --git a/src/mesh/group_manager.hh b/src/mesh/group_manager.hh index fd98c46b1..fb542763d 100644 --- a/src/mesh/group_manager.hh +++ b/src/mesh/group_manager.hh @@ -1,326 +1,356 @@ /** * @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)); } #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 and the associated node group - ElementGroup & createElementGroup(const std::string & group_name, - UInt dimension = _all_dimensions, - bool replace_group = false); - /// 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 &); 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/mesh/group_manager_inline_impl.cc b/src/mesh/group_manager_inline_impl.cc index bfb9d0f6a..418d6e627 100644 --- a/src/mesh/group_manager_inline_impl.cc +++ b/src/mesh/group_manager_inline_impl.cc @@ -1,205 +1,190 @@ /** * @file group_manager_inline_impl.cc * * @author Guillaume Anciaux * @author Dana Christen * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Sun Dec 03 2017 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "dumper_field.hh" #include "element_group.hh" #include "element_type_map_filter.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_nodal_field.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template class dump_type> dumper::Field * GroupManager::createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const ElementTypeMapArray * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") return this->createElementalField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); else return this->createElementalFilteredField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); } /* -------------------------------------------------------------------------- */ template class T2, template class, bool> class dump_type> dumper::Field * GroupManager::createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const ElementTypeMapArray * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") return this->createElementalField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); else return this->createElementalFilteredField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); } /* -------------------------------------------------------------------------- */ template class dump_type> ///< type of InternalMaterialField dumper::Field * GroupManager::createElementalField(const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const ElementTypeMapArray * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") return this->createElementalField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); else return this->createElementalFilteredField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); } /* -------------------------------------------------------------------------- */ template dumper::Field * GroupManager::createElementalField( const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, const ElementTypeMap & nb_data_per_elem) { const field_type * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name != "all") throw; dumper::Field * dumper = new dump_type(field, spatial_dimension, _not_ghost, kind); dumper->setNbDataPerElem(nb_data_per_elem); return dumper; } /* -------------------------------------------------------------------------- */ template dumper::Field * GroupManager::createElementalFilteredField( const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const field_type * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") throw; using T = typename field_type::type; ElementGroup & group = this->getElementGroup(group_name); UInt dim = group.getDimension(); if (dim != spatial_dimension) throw; const ElementTypeMapArray & elemental_filter = group.getElements(); auto * filtered = new ElementTypeMapArrayFilter(field, elemental_filter, nb_data_per_elem); dumper::Field * dumper = new dump_type(*filtered, dim, _not_ghost, kind); dumper->setNbDataPerElem(nb_data_per_elem); return dumper; } /* -------------------------------------------------------------------------- */ template class ftype> dumper::Field * GroupManager::createNodalField(const ftype * field, const std::string & group_name, UInt padding_size) { - if (field == nullptr) - return nullptr; - if (group_name == "all") { - using DumpType = typename dumper::NodalField; - auto * dumper = new DumpType(*field, 0, 0, nullptr); - dumper->setPadding(padding_size); - return dumper; - } else { - ElementGroup & group = this->getElementGroup(group_name); - const Array * nodal_filter = &(group.getNodes()); - using DumpType = typename dumper::NodalField; - auto * dumper = new DumpType(*field, 0, 0, nodal_filter); - dumper->setPadding(padding_size); - return dumper; - } - return nullptr; + return createStridedNodalField(field, group_name, 0, 0, padding_size); } /* -------------------------------------------------------------------------- */ template class ftype> dumper::Field * GroupManager::createStridedNodalField(const ftype * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size) { if (field == NULL) return nullptr; if (group_name == "all") { using DumpType = typename dumper::NodalField; auto * dumper = new DumpType(*field, size, stride, NULL); dumper->setPadding(padding_size); return dumper; } else { ElementGroup & group = this->getElementGroup(group_name); - const Array * nodal_filter = &(group.getNodes()); + const Array * nodal_filter = &(group.getNodeGroup().getNodes()); using DumpType = typename dumper::NodalField; auto * dumper = new DumpType(*field, size, stride, nodal_filter); dumper->setPadding(padding_size); return dumper; } return nullptr; } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif diff --git a/src/mesh/node_group.hh b/src/mesh/node_group.hh index f3b36be81..e991a1cf0 100644 --- a/src/mesh/node_group.hh +++ b/src/mesh/node_group.hh @@ -1,131 +1,132 @@ /** * @file node_group.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Nov 08 2017 * * @brief Node group definition * * @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_common.hh" #include "aka_memory.hh" #include "dumpable.hh" #include "mesh_filter.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_NODE_GROUP_HH__ #define __AKANTU_NODE_GROUP_HH__ namespace akantu { class NodeGroup : public Memory, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NodeGroup(const std::string & name, const Mesh & mesh, const std::string & id = "node_group", const MemoryID & memory_id = 0); ~NodeGroup() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: using const_node_iterator = Array::const_iterator; /// empty the node group void empty(); /// iterator to the beginning of the node group inline const_node_iterator begin() const; /// iterator to the end of the node group inline const_node_iterator end() const; /// add a node and give the local position through an iterator inline const_node_iterator add(UInt node, bool check_for_duplicate = true); /// remove a node inline void remove(UInt node); #ifndef SWIG inline decltype(auto) find(UInt node) const { return node_group.find(node); } #endif /// remove duplicated nodes void optimize(); /// append a group to current one void append(const NodeGroup & other_group); /// apply a filter on current node group template void applyNodeFilter(T & filter); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO_NOT_CONST(Nodes, node_group, Array &); AKANTU_GET_MACRO(Nodes, node_group, const Array &); AKANTU_GET_MACRO(Name, name, const std::string &); /// give the number of nodes in the current group inline UInt size() const; - UInt * storage() { return node_group.storage(); }; + //UInt * storage() { return node_group.storage(); }; + friend class GroupManager; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// name of the group std::string name; /// list of nodes in the group Array & node_group; /// reference to the mesh in question // const Mesh & mesh; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NodeGroup & _this) { _this.printself(stream); return stream; } } // akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "node_group_inline_impl.cc" #endif /* __AKANTU_NODE_GROUP_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh index dadf9e27e..2016ac6b0 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh @@ -1,801 +1,779 @@ /** * @file material_reinforcement_tmpl.hh * * @author Lucas Frerot * * @date creation: Wed Mar 25 2015 * @date last modification: Tue Feb 20 2018 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 2015-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_common.hh" #include "aka_voigthelper.hh" #include "material_reinforcement.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialReinforcement::MaterialReinforcement( EmbeddedInterfaceModel & model, const ID & id) : Mat(model, 1, model.getInterfaceMesh(), model.getFEEngine("EmbeddedInterfaceFEEngine"), id), emodel(model), gradu_embedded("gradu_embedded", *this, 1, model.getFEEngine("EmbeddedInterfaceFEEngine"), this->element_filter), directing_cosines("directing_cosines", *this, 1, model.getFEEngine("EmbeddedInterfaceFEEngine"), this->element_filter), pre_stress("pre_stress", *this, 1, model.getFEEngine("EmbeddedInterfaceFEEngine"), this->element_filter), area(1.0), shape_derivatives() { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initialize() { AKANTU_DEBUG_IN(); this->registerParam("area", area, _pat_parsable | _pat_modifiable, "Reinforcement cross-sectional area"); this->registerParam("pre_stress", pre_stress, _pat_parsable | _pat_modifiable, "Uniform pre-stress"); // this->unregisterInternal(this->stress); // Fool the AvgHomogenizingFunctor // stress.initialize(dim * dim); // Reallocate the element filter this->element_filter.initialize(this->emodel.getInterfaceMesh(), _spatial_dimension = 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialReinforcement::~MaterialReinforcement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initMaterial() { Mat::initMaterial(); gradu_embedded.initialize(dim * dim); pre_stress.initialize(1); /// We initialise the stuff that is not going to change during the simulation this->initFilters(); this->allocBackgroundShapeDerivatives(); this->initBackgroundShapeDerivatives(); this->initDirectingCosines(); } /* -------------------------------------------------------------------------- */ /// Initialize the filter for background elements template void MaterialReinforcement::initFilters() { for (auto gt : ghost_types) { for (auto && type : emodel.getInterfaceMesh().elementTypes(1, gt)) { std::string shaped_id = "filter"; if (gt == _ghost) shaped_id += ":ghost"; auto & background = background_filter(std::make_unique>( "bg_" + shaped_id, this->name), type, gt); auto & foreground = foreground_filter( std::make_unique>(shaped_id, this->name), type, gt); foreground->initialize(emodel.getMesh(), _nb_component = 1, _ghost_type = gt); background->initialize(emodel.getMesh(), _nb_component = 1, _ghost_type = gt); // Computing filters for (auto && bg_type : background->elementTypes(dim, gt)) { filterInterfaceBackgroundElements( (*foreground)(bg_type), (*background)(bg_type), bg_type, type, gt); } } } } /* -------------------------------------------------------------------------- */ /// Construct a filter for a (interface_type, background_type) pair template void MaterialReinforcement::filterInterfaceBackgroundElements( Array & foreground, Array & background, const ElementType & type, const ElementType & interface_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); foreground.resize(0); background.resize(0); Array & elements = emodel.getInterfaceAssociatedElements(interface_type, ghost_type); Array & elem_filter = this->element_filter(interface_type, ghost_type); for (auto & elem_id : elem_filter) { Element & elem = elements(elem_id); if (elem.type == type) { background.push_back(elem.element); foreground.push_back(elem_id); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ namespace detail { class BackgroundShapeDInitializer : public ElementTypeMapArrayInitializer { public: BackgroundShapeDInitializer(UInt spatial_dimension, FEEngine & engine, const ElementType & foreground_type, const ElementTypeMapArray & filter, const GhostType & ghost_type) : ElementTypeMapArrayInitializer( [](const ElementType & bgtype, const GhostType &) { return ShapeFunctions::getShapeDerivativesSize(bgtype); }, spatial_dimension, ghost_type, _ek_regular) { auto nb_quad = engine.getNbIntegrationPoints(foreground_type); // Counting how many background elements are affected by elements of // interface_type for (auto type : filter.elementTypes(this->spatial_dimension)) { // Inserting size array_size_per_bg_type(filter(type).size() * nb_quad, type, this->ghost_type); } } auto elementTypes() const -> decltype(auto) { return array_size_per_bg_type.elementTypes(); } UInt size(const ElementType & bgtype) const { return array_size_per_bg_type(bgtype, this->ghost_type); } protected: ElementTypeMap array_size_per_bg_type; }; -} +} // namespace detail /** * Background shape derivatives need to be stored per background element * types but also per embedded element type, which is why they are stored * in an ElementTypeMap *>. The outer ElementTypeMap * refers to the embedded types, and the inner refers to the background types. */ template void MaterialReinforcement::allocBackgroundShapeDerivatives() { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { for (auto && type : emodel.getInterfaceMesh().elementTypes(1, gt)) { std::string shaped_id = "embedded_shape_derivatives"; if (gt == _ghost) shaped_id += ":ghost"; auto & shaped_etma = shape_derivatives( std::make_unique>(shaped_id, this->name), type, gt); shaped_etma->initialize( detail::BackgroundShapeDInitializer( emodel.getSpatialDimension(), emodel.getFEEngine("EmbeddedInterfaceFEEngine"), type, *background_filter(type, gt), gt), 0, true); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initBackgroundShapeDerivatives() { AKANTU_DEBUG_IN(); for (auto interface_type : this->element_filter.elementTypes(this->spatial_dimension)) { for (auto type : background_filter(interface_type)->elementTypes(dim)) { computeBackgroundShapeDerivatives(interface_type, type, _not_ghost, this->element_filter(interface_type)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeBackgroundShapeDerivatives( const ElementType & interface_type, const ElementType & bg_type, GhostType ghost_type, const Array & filter) { auto & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); auto & engine = emodel.getFEEngine(); auto & interface_mesh = emodel.getInterfaceMesh(); const auto nb_nodes_elem_bg = Mesh::getNbNodesPerElement(bg_type); // const auto nb_strss = VoigtHelper::size; const auto nb_quads_per_elem = interface_engine.getNbIntegrationPoints(interface_type); Array quad_pos(0, dim, "interface_quad_pos"); interface_engine.interpolateOnIntegrationPoints(interface_mesh.getNodes(), quad_pos, dim, interface_type, ghost_type, filter); auto & background_shapesd = (*shape_derivatives(interface_type, ghost_type))(bg_type, ghost_type); auto & background_elements = (*background_filter(interface_type, ghost_type))(bg_type, ghost_type); auto & foreground_elements = (*foreground_filter(interface_type, ghost_type))(bg_type, ghost_type); auto shapesd_begin = background_shapesd.begin(dim, nb_nodes_elem_bg, nb_quads_per_elem); auto quad_begin = quad_pos.begin(dim, nb_quads_per_elem); for (auto && tuple : zip(background_elements, foreground_elements)) { UInt bg = std::get<0>(tuple), fg = std::get<1>(tuple); for (UInt i = 0; i < nb_quads_per_elem; ++i) { Matrix shapesd = Tensor3(shapesd_begin[fg])(i); Vector quads = Matrix(quad_begin[fg])(i); engine.computeShapeDerivatives(quads, bg, bg_type, shapesd, ghost_type); } } } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initDirectingCosines() { AKANTU_DEBUG_IN(); Mesh & mesh = emodel.getInterfaceMesh(); - Mesh::type_iterator type_it = mesh.firstType(1, _not_ghost); - Mesh::type_iterator type_end = mesh.lastType(1, _not_ghost); - const UInt voigt_size = VoigtHelper::size; directing_cosines.initialize(voigt_size); - for (; type_it != type_end; ++type_it) { - computeDirectingCosines(*type_it, _not_ghost); + for (auto && type : mesh.elementTypes(1, _not_ghost)) { + computeDirectingCosines(type, _not_ghost); // computeDirectingCosines(*type_it, _ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleStiffnessMatrix( GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = emodel.getInterfaceMesh(); - Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); - Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); - - for (; type_it != type_end; ++type_it) { - assembleStiffnessMatrix(*type_it, ghost_type); + for (auto && type : interface_mesh.elementTypes(1, _not_ghost)) { + assembleStiffnessMatrix(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleInternalForces( GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = emodel.getInterfaceMesh(); - Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); - Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); - - for (; type_it != type_end; ++type_it) { - this->assembleInternalForces(*type_it, ghost_type); + for (auto && type : interface_mesh.elementTypes(1, _not_ghost)) { + this->assembleInternalForces(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); - Mesh::type_iterator it = emodel.getInterfaceMesh().firstType(); - Mesh::type_iterator last_type = emodel.getInterfaceMesh().lastType(); - - for (; it != last_type; ++it) { - computeGradU(*it, ghost_type); - this->computeStress(*it, ghost_type); - addPrestress(*it, ghost_type); + Mesh & interface_mesh = emodel.getInterfaceMesh(); + for (auto && type : interface_mesh.elementTypes(_ghost_type = ghost_type)) { + computeGradU(type, ghost_type); + this->computeStress(type, ghost_type); + addPrestress(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::addPrestress(const ElementType & type, GhostType ghost_type) { auto & stress = this->stress(type, ghost_type); auto & sigma_p = this->pre_stress(type, ghost_type); for (auto && tuple : zip(stress, sigma_p)) { std::get<0>(tuple) += std::get<1>(tuple); } } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleInternalForces( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & mesh = emodel.getMesh(); - Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); - Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); - - for (; type_it != type_end; ++type_it) { - assembleInternalForcesInterface(type, *type_it, ghost_type); + for (auto && mesh_type : mesh.elementTypes(dim, ghost_type)) { + assembleInternalForcesInterface(type, mesh_type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Computes and assemble the residual. Residual in reinforcement is computed as: * * \f[ * \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s} * \f] */ template void MaterialReinforcement::assembleInternalForcesInterface( const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt voigt_size = VoigtHelper::size; FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); Array & elem_filter = this->element_filter(interface_type, ghost_type); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type); UInt nb_element = elem_filter.size(); UInt back_dof = dim * nodes_per_background_e; Array & shapesd = (*shape_derivatives(interface_type, ghost_type))( background_type, ghost_type); Array integrant(nb_quadrature_points * nb_element, back_dof, "integrant"); auto integrant_it = integrant.begin(back_dof); auto integrant_end = integrant.end(back_dof); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); auto C_it = directing_cosines(interface_type, ghost_type).begin(voigt_size); auto sigma_it = this->stress(interface_type, ghost_type).begin(); Matrix Bvoigt(voigt_size, back_dof); for (; integrant_it != integrant_end; ++integrant_it, ++B_it, ++C_it, ++sigma_it) { VoigtHelper::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt, nodes_per_background_e); Vector & C = *C_it; Vector & BtCt_sigma = *integrant_it; BtCt_sigma.mul(Bvoigt, C); BtCt_sigma *= *sigma_it * area; } Array residual_interface(nb_element, back_dof, "residual_interface"); interface_engine.integrate(integrant, residual_interface, back_dof, interface_type, ghost_type, elem_filter); integrant.resize(0); Array background_filter(nb_element, 1, "background_filter"); auto & filter = getBackgroundFilter(interface_type, background_type, ghost_type); emodel.getDOFManager().assembleElementalArrayLocalArray( residual_interface, emodel.getInternalForce(), background_type, ghost_type, -1., filter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeDirectingCosines( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = emodel.getInterfaceMesh(); const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const UInt steel_dof = dim * nb_nodes_per_element; const UInt voigt_size = VoigtHelper::size; const UInt nb_quad_points = emodel.getFEEngine("EmbeddedInterfaceFEEngine") .getNbIntegrationPoints(type, ghost_type); Array node_coordinates(this->element_filter(type, ghost_type).size(), steel_dof); this->emodel.getFEEngine().template extractNodalToElementField( interface_mesh, interface_mesh.getNodes(), node_coordinates, type, ghost_type, this->element_filter(type, ghost_type)); Array::matrix_iterator directing_cosines_it = directing_cosines(type, ghost_type).begin(1, voigt_size); Array::matrix_iterator node_coordinates_it = node_coordinates.begin(dim, nb_nodes_per_element); Array::matrix_iterator node_coordinates_end = node_coordinates.end(dim, nb_nodes_per_element); for (; node_coordinates_it != node_coordinates_end; ++node_coordinates_it) { for (UInt i = 0; i < nb_quad_points; i++, ++directing_cosines_it) { Matrix & nodes = *node_coordinates_it; Matrix & cosines = *directing_cosines_it; computeDirectingCosinesOnQuad(nodes, cosines); } } // Mauro: the directing_cosines internal is defined on the quadrature points // of each element AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleStiffnessMatrix( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & mesh = emodel.getMesh(); - Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); - Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); - - for (; type_it != type_end; ++type_it) { - assembleStiffnessMatrixInterface(type, *type_it, ghost_type); + for (auto && mesh_type : mesh.elementTypes(dim, ghost_type)) { + assembleStiffnessMatrixInterface(type, mesh_type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001) * \f[ * \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T * \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}} * \f] */ template void MaterialReinforcement::assembleStiffnessMatrixInterface( const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt voigt_size = VoigtHelper::size; FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); Array & elem_filter = this->element_filter(interface_type, ghost_type); Array & grad_u = gradu_embedded(interface_type, ghost_type); UInt nb_element = elem_filter.size(); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type); UInt back_dof = dim * nodes_per_background_e; UInt integrant_size = back_dof; grad_u.resize(nb_quadrature_points * nb_element); Array tangent_moduli(nb_element * nb_quadrature_points, 1, "interface_tangent_moduli"); this->computeTangentModuli(interface_type, tangent_moduli, ghost_type); Array & shapesd = (*shape_derivatives(interface_type, ghost_type))( background_type, ghost_type); Array integrant(nb_element * nb_quadrature_points, integrant_size * integrant_size, "B^t*C^t*D*C*B"); /// Temporary matrices for integrant product Matrix Bvoigt(voigt_size, back_dof); Matrix DCB(1, back_dof); Matrix CtDCB(voigt_size, back_dof); Array::scalar_iterator D_it = tangent_moduli.begin(); Array::scalar_iterator D_end = tangent_moduli.end(); Array::matrix_iterator C_it = directing_cosines(interface_type, ghost_type).begin(1, voigt_size); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator integrant_it = integrant.begin(integrant_size, integrant_size); for (; D_it != D_end; ++D_it, ++C_it, ++B_it, ++integrant_it) { Real & D = *D_it; Matrix & C = *C_it; Matrix & B = *B_it; Matrix & BtCtDCB = *integrant_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix(B, Bvoigt, nodes_per_background_e); DCB.mul(C, Bvoigt); DCB *= D * area; CtDCB.mul(C, DCB); BtCtDCB.mul(Bvoigt, CtDCB); } tangent_moduli.resize(0); Array K_interface(nb_element, integrant_size * integrant_size, "K_interface"); interface_engine.integrate(integrant, K_interface, integrant_size * integrant_size, interface_type, ghost_type, elem_filter); integrant.resize(0); // Mauro: Here K_interface contains the local stiffness matrices, // directing_cosines contains the information about the orientation // of the reinforcements, any rotation of the local stiffness matrix // can be done here auto & filter = getBackgroundFilter(interface_type, background_type, ghost_type); emodel.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", K_interface, background_type, ghost_type, _symmetric, filter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real MaterialReinforcement::getEnergy(const std::string & id) { AKANTU_DEBUG_IN(); if (id == "potential") { Real epot = 0.; this->computePotentialEnergyByElements(); - Mesh::type_iterator it = this->element_filter.firstType( - this->spatial_dimension), - end = this->element_filter.lastType( - this->spatial_dimension); - - for (; it != end; ++it) { + for (auto && type : this->element_filter.elementTypes(this->spatial_dimension)) { FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); epot += interface_engine.integrate( - this->potential_energy(*it, _not_ghost), *it, _not_ghost, - this->element_filter(*it, _not_ghost)); + this->potential_energy(type, _not_ghost), type, _not_ghost, + this->element_filter(type, _not_ghost)); epot *= area; } return epot; } AKANTU_DEBUG_OUT(); return 0; } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeGradU( const ElementType & interface_type, GhostType ghost_type) { // Looping over background types for (auto && bg_type : background_filter(interface_type, ghost_type)->elementTypes(dim)) { const UInt nodes_per_background_e = Mesh::getNbNodesPerElement(bg_type); const UInt voigt_size = VoigtHelper::size; auto & bg_shapesd = (*shape_derivatives(interface_type, ghost_type))(bg_type, ghost_type); auto & filter = getBackgroundFilter(interface_type, bg_type, ghost_type); Array disp_per_element(0, dim * nodes_per_background_e, "disp_elem"); FEEngine::extractNodalToElementField( emodel.getMesh(), emodel.getDisplacement(), disp_per_element, bg_type, ghost_type, filter); Matrix concrete_du(dim, dim); Matrix epsilon(dim, dim); Vector evoigt(voigt_size); for (auto && tuple : zip(make_view(disp_per_element, dim, nodes_per_background_e), make_view(bg_shapesd, dim, nodes_per_background_e), this->gradu(interface_type, ghost_type), make_view(this->directing_cosines(interface_type, ghost_type), voigt_size))) { auto & u = std::get<0>(tuple); auto & B = std::get<1>(tuple); auto & du = std::get<2>(tuple); auto & C = std::get<3>(tuple); concrete_du.mul(u, B); auto epsilon = 0.5 * (concrete_du + concrete_du.transpose()); strainTensorToVoigtVector(epsilon, evoigt); du = C.dot(evoigt); } } } /* -------------------------------------------------------------------------- */ /** * The structure of the directing cosines matrix is : * \f{eqnarray*}{ * C_{1,\cdot} & = & (l^2, m^2, n^2, mn, ln, lm) \\ * C_{i,j} & = & 0 * \f} * * with : * \f[ * (l, m, n) = \frac{1}{\|\frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}\|} \cdot * \frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s} * \f] */ template inline void MaterialReinforcement::computeDirectingCosinesOnQuad( const Matrix & nodes, Matrix & cosines) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(nodes.cols() == 2, "Higher order reinforcement elements not implemented"); const Vector a = nodes(0), b = nodes(1); Vector delta = b - a; Real sq_length = 0.; for (UInt i = 0; i < dim; i++) { sq_length += delta(i) * delta(i); } if (dim == 2) { cosines(0, 0) = delta(0) * delta(0); // l^2 cosines(0, 1) = delta(1) * delta(1); // m^2 cosines(0, 2) = delta(0) * delta(1); // lm } else if (dim == 3) { cosines(0, 0) = delta(0) * delta(0); // l^2 cosines(0, 1) = delta(1) * delta(1); // m^2 cosines(0, 2) = delta(2) * delta(2); // n^2 cosines(0, 3) = delta(1) * delta(2); // mn cosines(0, 4) = delta(0) * delta(2); // ln cosines(0, 5) = delta(0) * delta(1); // lm } cosines /= sq_length; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialReinforcement::stressTensorToVoigtVector( const Matrix & tensor, Vector & vector) { AKANTU_DEBUG_IN(); for (UInt i = 0; i < dim; i++) { vector(i) = tensor(i, i); } if (dim == 2) { vector(2) = tensor(0, 1); } else if (dim == 3) { vector(3) = tensor(1, 2); vector(4) = tensor(0, 2); vector(5) = tensor(0, 1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialReinforcement::strainTensorToVoigtVector( const Matrix & tensor, Vector & vector) { AKANTU_DEBUG_IN(); for (UInt i = 0; i < dim; i++) { vector(i) = tensor(i, i); } if (dim == 2) { vector(2) = 2 * tensor(0, 1); } else if (dim == 3) { vector(3) = 2 * tensor(1, 2); vector(4) = 2 * tensor(0, 2); vector(5) = 2 * tensor(0, 1); } AKANTU_DEBUG_OUT(); } -} // akantu +} // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.cc b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.cc index 29bd66b81..49caf656f 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_intersector.cc @@ -1,179 +1,174 @@ /** * @file embedded_interface_intersector.cc * * @author Lucas Frerot * * @date creation: Fri May 01 2015 * @date last modification: Tue Feb 20 2018 * * @brief Class that loads the interface from mesh and computes intersections * * @section LICENSE * * Copyright (©) 2015-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 "embedded_interface_intersector.hh" #include "mesh_segment_intersector.hh" /// Helper macro for types in the mesh. Creates an intersector and computes /// intersection queries #define INTERFACE_INTERSECTOR_CASE(dim, type) \ do { \ MeshSegmentIntersector intersector(this->mesh, interface_mesh); \ name_to_primitives_it = name_to_primitives_map.begin(); \ for (; name_to_primitives_it != name_to_primitives_end; \ ++name_to_primitives_it) { \ intersector.setPhysicalName(name_to_primitives_it->first); \ intersector.buildResultFromQueryList(name_to_primitives_it->second); \ } \ } while (0) #define INTERFACE_INTERSECTOR_CASE_2D(type) INTERFACE_INTERSECTOR_CASE(2, type) #define INTERFACE_INTERSECTOR_CASE_3D(type) INTERFACE_INTERSECTOR_CASE(3, type) namespace akantu { EmbeddedInterfaceIntersector::EmbeddedInterfaceIntersector( Mesh & mesh, const Mesh & primitive_mesh) : MeshGeomAbstract(mesh), interface_mesh(mesh.getSpatialDimension(), "interface_mesh"), primitive_mesh(primitive_mesh) { // Initiating mesh connectivity and data interface_mesh.addConnectivityType(_segment_2, _not_ghost); interface_mesh.addConnectivityType(_segment_2, _ghost); interface_mesh.registerElementalData("associated_element") .alloc(0, 1, _segment_2); interface_mesh.registerElementalData("physical_names") .alloc(0, 1, _segment_2); } EmbeddedInterfaceIntersector::~EmbeddedInterfaceIntersector() {} void EmbeddedInterfaceIntersector::constructData(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); const UInt dim = this->mesh.getSpatialDimension(); if (dim == 1) AKANTU_ERROR( "No embedded model in 1D. Deactivate intersection initialization"); Array * physical_names = NULL; try { physical_names = &const_cast &>( this->primitive_mesh.getData("physical_names", _segment_2)); } catch (debug::Exception & e) { AKANTU_ERROR("You must define physical names to reinforcements in " "order to use the embedded model"); throw e; } const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_segment_2); Array::const_vector_iterator connectivity = primitive_mesh.getConnectivity(_segment_2).begin(nb_nodes_per_element); Array::scalar_iterator names_it = physical_names->begin(), names_end = physical_names->end(); std::map> name_to_primitives_map; // Loop over the physical names and register segment lists in // name_to_primitives_map for (; names_it != names_end; ++names_it) { UInt element_id = names_it - physical_names->begin(); const Vector el_connectivity = connectivity[element_id]; K::Segment_3 segment = this->createSegment(el_connectivity); name_to_primitives_map[*names_it].push_back(segment); } // Loop over the background types of the mesh - Mesh::type_iterator type_it = this->mesh.firstType(dim, _not_ghost), - type_end = this->mesh.lastType(dim, _not_ghost); - std::map>::iterator name_to_primitives_it, name_to_primitives_end = name_to_primitives_map.end(); - for (; type_it != type_end; ++type_it) { + for (auto type : this->mesh.elementTypes(dim, _not_ghost)) { // Used in AKANTU_BOOST_ELEMENT_SWITCH - ElementType type = *type_it; - AKANTU_DEBUG_INFO("Computing intersections with background element type " << type); switch (dim) { case 1: break; case 2: // Compute intersections for supported 2D elements AKANTU_BOOST_ELEMENT_SWITCH(INTERFACE_INTERSECTOR_CASE_2D, (_triangle_3)(_triangle_6)); break; case 3: // Compute intersections for supported 3D elements AKANTU_BOOST_ELEMENT_SWITCH(INTERFACE_INTERSECTOR_CASE_3D, (_tetrahedron_4)); break; } } AKANTU_DEBUG_OUT(); } K::Segment_3 EmbeddedInterfaceIntersector::createSegment(const Vector & connectivity) { AKANTU_DEBUG_IN(); K::Point_3 *source = NULL, *target = NULL; const Array & nodes = this->primitive_mesh.getNodes(); if (this->mesh.getSpatialDimension() == 2) { source = new K::Point_3(nodes(connectivity(0), 0), nodes(connectivity(0), 1), 0.); target = new K::Point_3(nodes(connectivity(1), 0), nodes(connectivity(1), 1), 0.); } else if (this->mesh.getSpatialDimension() == 3) { source = new K::Point_3(nodes(connectivity(0), 0), nodes(connectivity(0), 1), nodes(connectivity(0), 2)); target = new K::Point_3(nodes(connectivity(1), 0), nodes(connectivity(1), 1), nodes(connectivity(1), 2)); } K::Segment_3 segment(*source, *target); delete source; delete target; AKANTU_DEBUG_OUT(); return segment; } } // namespace akantu #undef INTERFACE_INTERSECTOR_CASE #undef INTERFACE_INTERSECTOR_CASE_2D #undef INTERFACE_INTERSECTOR_CASE_3D diff --git a/test/test_io/test_dumper/test_dumper.cc b/test/test_io/test_dumper/test_dumper.cc index 93a39389e..e4b5ccc35 100644 --- a/test/test_io/test_dumper/test_dumper.cc +++ b/test/test_io/test_dumper/test_dumper.cc @@ -1,156 +1,156 @@ /** * @file test_dumper.cc * * @author David Simon Kammer * * @date creation: Tue Sep 02 2014 * @date last modification: Mon Jan 22 2018 * * @brief test dumper * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" #include "dumper_variable.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize("input_file.dat", argc, argv); UInt spatial_dimension = 3; Mesh mesh(spatial_dimension); mesh.read("test_dumper.msh"); SolidMechanicsModel model(mesh); auto && mat_selector = std::make_shared>("physical_names", model); model.setMaterialSelector(mat_selector); model.initFull(); model.assembleInternalForces(); Real time_step = 0.1; const Array & coord = mesh.getNodes(); Array & disp = model.getDisplacement(); Array & bound = model.getBlockedDOFs(); for (UInt n = 0; n < mesh.getNbNodes(); ++n) { Real dist = 0.; for (UInt d = 0; d < spatial_dimension; ++d) { dist += coord(n, d) * coord(n, d); } dist = sqrt(dist); for (UInt d = 0; d < spatial_dimension; ++d) { disp(n, d) = (d + 1) * dist; bound(n, d) = bool((n % 2) + d); } } // dump boundary bottom as reference model.setGroupDirectory("paraview", "Bottom"); model.setGroupBaseName("paraview_bottom", "Bottom"); model.addDumpGroupField("displacement", "Bottom"); model.addDumpGroupField("blocked_dofs", "Bottom"); UInt nbp = 3; DumperParaview prvdumper("paraview_bottom_parallel", "paraview", false); iohelper::Dumper & prvdpr = prvdumper.getDumper(); for (UInt p = 0; p < nbp; ++p) { prvdpr.setParallelContext(p, nbp, 0); if (p != 0) { prvdumper.unRegisterField("connectivities"); prvdumper.unRegisterField("element_type"); prvdumper.unRegisterField("positions"); prvdumper.unRegisterField("displacement"); } prvdumper.registerFilteredMesh(mesh, mesh.getElementGroup("Bottom").getElements(), - mesh.getElementGroup("Bottom").getNodes()); + mesh.getElementGroup("Bottom").getNodeGroup().getNodes()); prvdumper.registerField("displacement", new dumper::NodalField( model.getDisplacement(), 0, 0, - &(mesh.getElementGroup("Bottom").getNodes()))); + &(mesh.getElementGroup("Bottom").getNodeGroup().getNodes()))); prvdumper.dump(0); } DumperText txtdumper("text_bottom", iohelper::_tdm_csv); txtdumper.setDirectory("paraview"); txtdumper.setPrecision(8); txtdumper.setTimeStep(time_step); txtdumper.registerFilteredMesh(mesh, mesh.getElementGroup("Bottom").getElements(), - mesh.getElementGroup("Bottom").getNodes()); + mesh.getElementGroup("Bottom").getNodeGroup().getNodes()); txtdumper.registerField("displacement", new dumper::NodalField( model.getDisplacement(), 0, 0, - &(mesh.getElementGroup("Bottom").getNodes()))); + &(mesh.getElementGroup("Bottom").getNodeGroup().getNodes()))); txtdumper.registerField("blocked_dofs", new dumper::NodalField( model.getBlockedDOFs(), 0, 0, - &(mesh.getElementGroup("Bottom").getNodes()))); + &(mesh.getElementGroup("Bottom").getNodeGroup().getNodes()))); Real pot_energy = 1.2345567891; Vector gforces(2, 1.); txtdumper.registerVariable("potential_energy", new dumper::Variable(pot_energy)); txtdumper.registerVariable("global_forces", new dumper::Variable>(gforces)); // dump a first time before the main loop model.dumpGroup("Bottom"); txtdumper.dump(); Real time = 0.; for (UInt i = 1; i < 5; ++i) { pot_energy += 2.; gforces(0) += 0.1; gforces(1) += 0.2; // pre -> cor // increment time after all steps of integration time += time_step; // dump after time increment if (i % 2 == 0) { txtdumper.dump(time, i); model.dumpGroup("Bottom"); // parallel test for (UInt p = 0; p < nbp; ++p) { prvdpr.setParallelContext(p, nbp, 0); prvdumper.dump(i); } } } finalize(); return EXIT_SUCCESS; } diff --git a/test/test_mesh_utils/test_mesh_iterators.cc b/test/test_mesh_utils/test_mesh_iterators.cc index bd27aa63b..2bf1ad4e9 100644 --- a/test/test_mesh_utils/test_mesh_iterators.cc +++ b/test/test_mesh_utils/test_mesh_iterators.cc @@ -1,70 +1,73 @@ /** * @file test_mesh_iterators.cc * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Aug 16 2017 * * @brief Test the mesh iterators * * @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_iterators.hh" #include "element_group.hh" #include "mesh.hh" #include "mesh_iterators.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize(argc, argv); Mesh mesh(3); const Mesh & cmesh = mesh; mesh.read("iterators_mesh.msh"); + std::cout << "ElementGroups" << std::endl; for (auto && element_group : ElementGroupsIterable(mesh)) { - std::cout << element_group.getName() << std::endl; + std::cout << element_group.getName() << " " << element_group.getDimension() << std::endl; } + std::cout << "NodeGroups" << std::endl; for (auto && node_group : NodeGroupsIterable(cmesh)) { std::cout << node_group.getName() << std::endl; } + std::cout << "enumerate(ElementGroups)" << std::endl; for (auto && element_group : enumerate(ElementGroupsIterable(mesh))) { std::cout << std::get<0>(element_group) << " " << std::get<1>(element_group).getName() << std::endl; } // for (auto && node_group : // counting(NodeGroupsIterable(cmesh))) { // std::cout << std::get<0>(node_group) << " " << // std::get<1>(node_group).getName() << std::endl; // } finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/patch_tests/CMakeLists.txt b/test/test_model/patch_tests/CMakeLists.txt index 781b4c8da..32bec6111 100644 --- a/test/test_model/patch_tests/CMakeLists.txt +++ b/test/test_model/patch_tests/CMakeLists.txt @@ -1,124 +1,124 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author David Simon Kammer # @author Nicolas Richart # # @date creation: Fri Oct 22 2010 # @date last modification: Thu Feb 08 2018 # # @brief configuration for patch tests # # @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 . # # @section DESCRIPTION # #=============================================================================== add_subdirectory(data) register_gtest_sources(SOURCES patch_test_linear_elastic_explicit.cc PACKAGE solid_mechanics FILES_TO_COPY data/material_check_stress_plane_strain.dat data/material_check_stress_plane_stress.dat) register_gtest_sources(SOURCES patch_test_linear_elastic_implicit.cc PACKAGE solid_mechanics implicit FILES_TO_COPY data/material_check_stress_plane_strain.dat data/material_check_stress_plane_stress.dat) register_gtest_sources(SOURCES patch_test_linear_anisotropic.cc PACKAGE solid_mechanics lapack FILES_TO_COPY data/material_anisotropic_1.dat data/material_anisotropic_2.dat data/material_anisotropic_3.dat ) register_gtest_sources(SOURCES test_lumped_mass.cc PACKAGE solid_mechanics FILES_TO_COPY data/material_lumped.dat) register_gtest_sources( SOURCES patch_test_linear_heat_transfer_explicit.cc FILES_TO_COPY data/heat_transfer_input.dat PACKAGE heat_transfer) register_gtest_sources( SOURCES patch_test_linear_heat_transfer_static.cc FILES_TO_COPY data/heat_transfer_input.dat PACKAGE heat_transfer implicit) register_gtest_sources( SOURCES patch_test_linear_heat_transfer_implicit.cc FILES_TO_COPY data/heat_transfer_input.dat PACKAGE heat_transfer implicit ) register_gtest_test(patch_test_linear FILES_TO_COPY ${PATCH_TEST_MESHES}) register_test(test_linear_elastic_explicit_python SCRIPT patch_test_linear_elastic_explicit.py PYTHON PACKAGE python_interface FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_solid_mechanics_fixture.py FILES_TO_COPY ${PATCH_TEST_MESHES}) register_test(test_linear_elastic_static_python SCRIPT patch_test_linear_elastic_static.py PYTHON PACKAGE python_interface implicit FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_solid_mechanics_fixture.py) register_test(test_linear_anisotropic_explicit_python SCRIPT patch_test_linear_anisotropic_explicit.py PYTHON UNSTABLE PACKAGE python_interface lapack FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_solid_mechanics_fixture.py - FILES_TO_COPY data/material_anisotropic.dat) + FILES_TO_COPY data/material_anisotropic_3.dat) register_test(patch_test_linear_heat_transfer_explicit_python SCRIPT patch_test_linear_heat_transfer_explicit.py PYTHON PACKAGE python_interface heat_transfer FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_heat_transfer_fixture.py FILES_TO_COPY data/heat_transfer_input.dat) register_test(patch_test_linear_heat_transfer_static_python SCRIPT patch_test_linear_heat_transfer_static.py PYTHON PACKAGE python_interface heat_transfer implicit FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_heat_transfer_fixture.py FILES_TO_COPY data/heat_transfer_input.dat) register_test(patch_test_linear_heat_transfer_implicit_python SCRIPT patch_test_linear_heat_transfer_implicit.py PYTHON PACKAGE python_interface heat_transfer implicit FILES_TO_COPY patch_test_linear_fixture.py FILES_TO_COPY patch_test_linear_heat_transfer_fixture.py FILES_TO_COPY data/heat_transfer_input.dat) diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc index 510aff23b..56e283ba9 100644 --- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc +++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc @@ -1,234 +1,234 @@ /** * @file test_embedded_interface_model_prestress.cc * * @author Lucas Frerot * * @date creation: Tue Apr 28 2015 * @date last modification: Tue Feb 20 2018 * * @brief Embedded model test for prestressing (bases on stress norm) * * @section LICENSE * * Copyright (©) 2015-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_common.hh" #include "embedded_interface_model.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; #define YG 0.483644859 #define I_eq 0.012488874 #define A_eq (1e-2 + 1. / 7. * 1.) /* -------------------------------------------------------------------------- */ struct StressSolution : public BC::Neumann::FromHigherDim { Real M; Real I; Real yg; Real pre_stress; StressSolution(UInt dim, Real M, Real I, Real yg = 0, Real pre_stress = 0) : BC::Neumann::FromHigherDim(Matrix(dim, dim)), M(M), I(I), yg(yg), pre_stress(pre_stress) {} virtual ~StressSolution() {} void operator()(const IntegrationPoint & /*quad_point*/, Vector & dual, const Vector & coord, const Vector & normals) const { UInt dim = coord.size(); if (dim < 2) AKANTU_ERROR("Solution not valid for 1D"); Matrix stress(dim, dim); stress.clear(); stress(0, 0) = this->stress(coord(1)); dual.mul(stress, normals); } inline Real stress(Real height) const { return -M / I * (height - yg) + pre_stress; } inline Real neutral_axis() const { return -I * pre_stress / M + yg; } }; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("prestress.dat", argc, argv); debug::setDebugLevel(dblError); Math::setTolerance(1e-6); const UInt dim = 2; /* -------------------------------------------------------------------------- */ Mesh mesh(dim); mesh.read("embedded_mesh_prestress.msh"); // mesh.createGroupsFromMeshData("physical_names"); Mesh reinforcement_mesh(dim, "reinforcement_mesh"); try { reinforcement_mesh.read("embedded_mesh_prestress_reinforcement.msh"); } catch (debug::Exception & e) { } // reinforcement_mesh.createGroupsFromMeshData("physical_names"); EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim); model.initFull(EmbeddedInterfaceModelOptions(_static)); /* -------------------------------------------------------------------------- */ /* Computation of analytical residual */ /* -------------------------------------------------------------------------- */ /* * q = 1000 N/m * L = 20 m * a = 1 m */ Real steel_area = model.getMaterial("reinforcement").get("area"); Real pre_stress = model.getMaterial("reinforcement").get("pre_stress"); Real stress_norm = 0.; StressSolution *concrete_stress = nullptr, *steel_stress = nullptr; Real pre_force = pre_stress * steel_area; Real pre_moment = -pre_force * (YG - 0.25); Real neutral_axis = YG - I_eq / A_eq * pre_force / pre_moment; concrete_stress = new StressSolution(dim, pre_moment, 7. * I_eq, YG, -pre_force / (7. * A_eq)); steel_stress = new StressSolution(dim, pre_moment, I_eq, YG, pre_stress - pre_force / A_eq); stress_norm = std::abs(concrete_stress->stress(1)) * (1 - neutral_axis) * 0.5 + std::abs(concrete_stress->stress(0)) * neutral_axis * 0.5 + std::abs(steel_stress->stress(0.25)) * steel_area; model.applyBC(*concrete_stress, "XBlocked"); auto end_node = *mesh.getElementGroup("EndNode").getNodeGroup().begin(); Vector end_node_force = model.getExternalForce().begin(dim)[end_node]; end_node_force(0) += steel_stress->stress(0.25) * steel_area; Array analytical_residual(mesh.getNbNodes(), dim, "analytical_residual"); analytical_residual.copy(model.getExternalForce()); model.getExternalForce().clear(); delete concrete_stress; delete steel_stress; /* -------------------------------------------------------------------------- */ model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "XBlocked"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "YBlocked"); try { model.solveStep(); - } catch (debug::Exception e) { + } catch (debug::Exception & e) { std::cerr << e.what() << std::endl; return EXIT_FAILURE; } /* -------------------------------------------------------------------------- */ /* Computation of FEM residual norm */ /* -------------------------------------------------------------------------- */ ElementGroup & xblocked = mesh.getElementGroup("XBlocked"); NodeGroup & boundary_nodes = xblocked.getNodeGroup(); NodeGroup::const_node_iterator nodes_it = boundary_nodes.begin(), nodes_end = boundary_nodes.end(); model.assembleInternalForces(); Array residual(mesh.getNbNodes(), dim, "my_residual"); residual.copy(model.getInternalForce()); residual -= model.getExternalForce(); auto com_res = residual.begin(dim); auto position = mesh.getNodes().begin(dim); Real res_sum = 0.; UInt lower_node = -1; UInt upper_node = -1; Real lower_dist = 1; Real upper_dist = 1; for (; nodes_it != nodes_end; ++nodes_it) { UInt node_number = *nodes_it; const Vector res = com_res[node_number]; const Vector pos = position[node_number]; if (!Math::are_float_equal(pos(1), 0.25)) { if ((std::abs(pos(1) - 0.25) < lower_dist) && (pos(1) < 0.25)) { lower_dist = std::abs(pos(1) - 0.25); lower_node = node_number; } if ((std::abs(pos(1) - 0.25) < upper_dist) && (pos(1) > 0.25)) { upper_dist = std::abs(pos(1) - 0.25); upper_node = node_number; } } for (UInt i = 0; i < dim; i++) { if (!Math::are_float_equal(pos(1), 0.25)) { res_sum += std::abs(res(i)); } } } const Vector upper_res = com_res[upper_node], lower_res = com_res[lower_node]; const Vector end_node_res = com_res[end_node]; Vector delta = upper_res - lower_res; delta *= lower_dist / (upper_dist + lower_dist); Vector concrete_residual = lower_res + delta; Vector steel_residual = end_node_res - concrete_residual; for (UInt i = 0; i < dim; i++) { res_sum += std::abs(concrete_residual(i)); res_sum += std::abs(steel_residual(i)); } Real relative_error = std::abs(res_sum - stress_norm) / stress_norm; if (relative_error > 1e-3) { std::cerr << "Relative error = " << relative_error << std::endl; return EXIT_FAILURE; } finalize(); return 0; } diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_elasto_plastic_linear_isotropic_hardening/test_material_elasto_plastic_linear_isotropic_hardening.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_elasto_plastic_linear_isotropic_hardening/test_material_elasto_plastic_linear_isotropic_hardening.cc index 90fc154a9..abad2e026 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_elasto_plastic_linear_isotropic_hardening/test_material_elasto_plastic_linear_isotropic_hardening.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_elasto_plastic_linear_isotropic_hardening/test_material_elasto_plastic_linear_isotropic_hardening.cc @@ -1,91 +1,91 @@ /** * @file test_material_elasto_plastic_linear_isotropic_hardening.cc * * @author Jaehyun Cho * @author Lucas Frerot * * @date creation: Sun Oct 19 2014 * @date last modification: Mon Sep 11 2017 * * @brief test for material type elasto plastic linear isotropic hardening * using tension-compression test * * @section LICENSE * * Copyright (©) 2015-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 "non_linear_solver.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("test_material_elasto_plastic_linear_isotropic_hardening.dat", argc, argv); const UInt spatial_dimension = 2; const Real u_increment = 0.1; const UInt steps = 20; Mesh mesh(spatial_dimension); mesh.read("test_material_elasto_plastic_linear_isotropic_hardening.msh"); SolidMechanicsModel model(mesh); model.initFull(_analysis_method = _static); auto & solver = model.getNonLinearSolver("static"); solver.set("max_iterations", 300); solver.set("threshold", 1e-5); model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "left"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "bottom"); std::cout.precision(4); for (UInt i = 0; i < steps; ++i) { model.applyBC(BC::Dirichlet::FixedValue(i * u_increment, _x), "right"); try { model.solveStep(); } catch (debug::NLSNotConvergedException & e) { std::cout << e.niter << " " << e.error << std::endl; throw; } Real strainxx = i * u_increment / 10.; - const Array & edge_nodes = mesh.getElementGroup("right").getNodes(); + const Array & edge_nodes = mesh.getElementGroup("right").getNodeGroup().getNodes(); Array & residual = model.getInternalForce(); Real reaction = 0; for (UInt n = 0; n < edge_nodes.size(); n++) { reaction -= residual(edge_nodes(n), 0); } std::cout << strainxx << "," << reaction << std::endl; } finalize(); return 0; }