diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc index 30b05d7fc..f90facf28 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc @@ -1,558 +1,555 @@ /** * @file ntn_base_contact.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of ntn base 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 . * */ /* -------------------------------------------------------------------------- */ #include "ntn_base_contact.hh" #include "dof_manager_default.hh" #include "dumpable_inline_impl.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" #include "element_synchronizer.hh" #include "mesh_utils.hh" #include "non_linear_solver_lumped.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ // NTNContactSynchElementFilter::NTNContactSynchElementFilter( // NTNBaseContact & contact) // : contact(contact), // connectivity(contact.getModel().getMesh().getConnectivities()) { // AKANTU_DEBUG_IN(); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ // bool NTNContactSynchElementFilter::operator()(const Element & e) { // AKANTU_DEBUG_IN(); // ElementType type = e.type; // UInt element = e.element; // GhostType ghost_type = e.ghost_type; // // loop over all nodes of this element // bool need_element = false; // UInt nb_nodes = Mesh::getNbNodesPerElement(type); // for (UInt n = 0; n < nb_nodes; ++n) { // UInt nn = this->connectivity(type, ghost_type)(element, n); // // if one nodes is in this contact, we need this element // if (this->contact.getNodeIndex(nn) >= 0) { // need_element = true; // break; // } // } // AKANTU_DEBUG_OUT(); // return need_element; // } /* -------------------------------------------------------------------------- */ NTNBaseContact::NTNBaseContact(SolidMechanicsModel & model, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), Dumpable(), model(model), slaves(0, 1, 0, id + ":slaves", std::numeric_limits::quiet_NaN(), "slaves"), normals(0, model.getSpatialDimension(), 0, id + ":normals", std::numeric_limits::quiet_NaN(), "normals"), contact_pressure( 0, model.getSpatialDimension(), 0, id + ":contact_pressure", std::numeric_limits::quiet_NaN(), "contact_pressure"), is_in_contact(0, 1, false, id + ":is_in_contact", false, "is_in_contact"), lumped_boundary_slaves(0, 1, 0, id + ":lumped_boundary_slaves", std::numeric_limits::quiet_NaN(), "lumped_boundary_slaves"), impedance(0, 1, 0, id + ":impedance", std::numeric_limits::quiet_NaN(), "impedance"), contact_surfaces(), slave_elements("slave_elements", id, memory_id), node_to_elements() { AKANTU_DEBUG_IN(); FEEngine & boundary_fem = this->model.getFEEngineBoundary(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { boundary_fem.initShapeFunctions(*gt); } Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = this->model.getSpatialDimension(); this->slave_elements.initialize(mesh, _spatial_dimension = spatial_dimension - 1); MeshUtils::buildNode2Elements(mesh, this->node_to_elements, spatial_dimension - 1); this->registerDumper("text_all", id, true); this->addDumpFilteredMesh(mesh, slave_elements, slaves.getArray(), spatial_dimension - 1, _not_ghost, _ek_regular); // parallelisation this->synch_registry = std::make_unique(); this->synch_registry->registerDataAccessor(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ NTNBaseContact::~NTNBaseContact() = default; /* -------------------------------------------------------------------------- */ void NTNBaseContact::initParallel() { AKANTU_DEBUG_IN(); this->synchronizer = std::make_unique( this->model.getMesh().getElementSynchronizer()); this->synchronizer->filterScheme([&](auto && element) { // loop over all nodes of this element Vector conn = const_cast(this->model.getMesh()) .getConnectivity(element); for (auto & node : conn) { // if one nodes is in this contact, we need this element if (this->getNodeIndex(node) >= 0) { return true; } } return false; }); this->synch_registry->registerSynchronizer(*(this->synchronizer), _gst_cf_nodal); this->synch_registry->registerSynchronizer(*(this->synchronizer), _gst_cf_incr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::findBoundaryElements( const Array & interface_nodes, ElementTypeMapArray & elements) { AKANTU_DEBUG_IN(); // add connected boundary elements that have all nodes on this contact for (const auto & node : interface_nodes) { for (const auto & element : this->node_to_elements.getRow(node)) { Vector conn = const_cast(this->model.getMesh()) .getConnectivity(element); auto nb_nodes = conn.size(); decltype(nb_nodes) nb_found_nodes = 0; for (auto & nn : conn) { if (interface_nodes.find(nn) != UInt(-1)) { nb_found_nodes++; } else { break; } } // this is an element between all contact nodes // and is not already in the elements if ((nb_found_nodes == nb_nodes) && (elements(element.type, element.ghost_type).find(element.element) == UInt(-1))) { elements(element.type, element.ghost_type).push_back(element.element); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::addSplitNode(UInt node) { AKANTU_DEBUG_IN(); UInt dim = this->model.getSpatialDimension(); // add to node arrays this->slaves.push_back(node); // set contact as false this->is_in_contact.push_back(false); // before initializing // set contact pressure, normal, lumped_boundary to Nan this->contact_pressure.push_back(std::numeric_limits::quiet_NaN()); this->impedance.push_back(std::numeric_limits::quiet_NaN()); this->lumped_boundary_slaves.push_back( std::numeric_limits::quiet_NaN()); Vector nan_normal(dim, std::numeric_limits::quiet_NaN()); this->normals.push_back(nan_normal); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::registerSynchronizedArray(SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->slaves.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); this->slaves.dumpRestartFile(file_name); this->normals.dumpRestartFile(file_name); this->is_in_contact.dumpRestartFile(file_name); this->contact_pressure.dumpRestartFile(file_name); this->lumped_boundary_slaves.dumpRestartFile(file_name); this->impedance.dumpRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); this->slaves.readRestartFile(file_name); this->normals.readRestartFile(file_name); this->is_in_contact.readRestartFile(file_name); this->contact_pressure.readRestartFile(file_name); this->lumped_boundary_slaves.readRestartFile(file_name); this->impedance.readRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt NTNBaseContact::getNbNodesInContact() const { AKANTU_DEBUG_IN(); UInt nb_contact = 0; UInt nb_nodes = this->getNbContactNodes(); const Mesh & mesh = this->model.getMesh(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(this->slaves(n)); bool is_pbc_slave_node = mesh.isPeriodicSlave(this->slaves(n)); if (is_local_node && !is_pbc_slave_node && this->is_in_contact(n)) { nb_contact++; } } mesh.getCommunicator().allReduce(nb_contact, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return nb_contact; } /* -------------------------------------------------------------------------- */ void NTNBaseContact::updateInternalData() { AKANTU_DEBUG_IN(); updateNormals(); updateLumpedBoundary(); updateImpedance(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::updateLumpedBoundary() { AKANTU_DEBUG_IN(); this->internalUpdateLumpedBoundary(this->slaves.getArray(), this->slave_elements, this->lumped_boundary_slaves); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::internalUpdateLumpedBoundary( const Array & nodes, const ElementTypeMapArray & elements, SynchronizedArray & boundary) { AKANTU_DEBUG_IN(); // set all values in lumped_boundary to zero boundary.clear(); UInt dim = this->model.getSpatialDimension(); // UInt nb_contact_nodes = getNbContactNodes(); const FEEngine & boundary_fem = this->model.getFEEngineBoundary(); const Mesh & mesh = this->model.getMesh(); - for (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - Mesh::type_iterator it = mesh.firstType(dim - 1, *gt); - Mesh::type_iterator last = mesh.lastType(dim - 1, *gt); - for (; it != last; ++it) { - UInt nb_elements = mesh.getNbElement(*it, *gt); - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); - const Array & connectivity = mesh.getConnectivity(*it, *gt); + for(auto ghost_type : ghost_types) { + for (auto & type : mesh.elementTypes(dim - 1, ghost_type)) { + UInt nb_elements = mesh.getNbElement(type, ghost_type); + UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); + const Array & connectivity = mesh.getConnectivity(type, ghost_type); // get shapes and compute integral - const Array & shapes = boundary_fem.getShapes(*it, *gt); + const Array & shapes = boundary_fem.getShapes(type, ghost_type); Array area(nb_elements, nb_nodes_per_element); - boundary_fem.integrate(shapes, area, nb_nodes_per_element, *it, *gt); + boundary_fem.integrate(shapes, area, nb_nodes_per_element, type, ghost_type); if (this->contact_surfaces.size() == 0) { AKANTU_DEBUG_WARNING( "No surfaces in ntn base contact." << " You have to define the lumped boundary by yourself."); } - Array::const_iterator elem_it = (elements)(*it, *gt).begin(); + Array::const_iterator elem_it = (elements)(type, ghost_type).begin(); Array::const_iterator elem_it_end = - (elements)(*it, *gt).end(); + (elements)(type, ghost_type).end(); // loop over contact nodes for (; elem_it != elem_it_end; ++elem_it) { for (UInt q = 0; q < nb_nodes_per_element; ++q) { UInt node = connectivity(*elem_it, q); UInt node_index = nodes.find(node); AKANTU_DEBUG_ASSERT(node_index != UInt(-1), "Could not find node " << node << " in the array!"); Real area_to_add = area(*elem_it, q); boundary(node_index) += area_to_add; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::computeContactPressure() { AKANTU_DEBUG_IN(); UInt dim = this->model.getSpatialDimension(); Real delta_t = this->model.getTimeStep(); UInt nb_contact_nodes = getNbContactNodes(); AKANTU_DEBUG_ASSERT(delta_t > 0., "Cannot compute contact pressure if no time step is set"); // synchronize data this->synch_registry->synchronize(_gst_cf_nodal); auto && dof_manager = dynamic_cast(model.getDOFManager()); const auto & b = dof_manager.getResidual(); Array acceleration(b.size(), dim); const auto & blocked_dofs = dof_manager.getGlobalBlockedDOFs(); const auto & A = dof_manager.getLumpedMatrix("M"); // pre-compute the acceleration // (not increment acceleration, because residual is still Kf) NonLinearSolverLumped::solveLumped(A, acceleration, b, blocked_dofs, this->model.getF_M2A()); // compute relative normal fields of displacement, velocity and acceleration Array r_disp(0, 1); Array r_velo(0, 1); Array r_acce(0, 1); Array r_old_acce(0, 1); computeNormalGap(r_disp); // computeRelativeNormalField(this->model.getCurrentPosition(), r_disp); computeRelativeNormalField(this->model.getVelocity(), r_velo); computeRelativeNormalField(acceleration, r_acce); computeRelativeNormalField(this->model.getAcceleration(), r_old_acce); AKANTU_DEBUG_ASSERT(r_disp.size() == nb_contact_nodes, "computeRelativeNormalField does not give back arrays " << "size == nb_contact_nodes. nb_contact_nodes = " << nb_contact_nodes << " | array size = " << r_disp.size()); // compute gap array for all nodes Array gap(nb_contact_nodes, 1); Real * gap_p = gap.storage(); Real * r_disp_p = r_disp.storage(); Real * r_velo_p = r_velo.storage(); Real * r_acce_p = r_acce.storage(); Real * r_old_acce_p = r_old_acce.storage(); for (UInt i = 0; i < nb_contact_nodes; ++i) { *gap_p = *r_disp_p + delta_t * *r_velo_p + delta_t * delta_t * *r_acce_p - 0.5 * delta_t * delta_t * *r_old_acce_p; // increment pointers gap_p++; r_disp_p++; r_velo_p++; r_acce_p++; r_old_acce_p++; } // check if gap is negative -> is in contact for (UInt n = 0; n < nb_contact_nodes; ++n) { if (gap(n) <= 0.) { for (UInt d = 0; d < dim; ++d) { this->contact_pressure(n, d) = this->impedance(n) * gap(n) / (2 * delta_t) * this->normals(n, d); } this->is_in_contact(n) = true; } else { for (UInt d = 0; d < dim; ++d) this->contact_pressure(n, d) = 0.; this->is_in_contact(n) = false; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::applyContactPressure() { AKANTU_DEBUG_IN(); UInt nb_contact_nodes = getNbContactNodes(); UInt dim = this->model.getSpatialDimension(); Array & residual = this->model.getInternalForce(); for (UInt n = 0; n < nb_contact_nodes; ++n) { UInt slave = this->slaves(n); for (UInt d = 0; d < dim; ++d) { // residual(master,d) += this->lumped_boundary(n,0) * // this->contact_pressure(n,d); residual(slave, d) -= this->lumped_boundary_slaves(n) * this->contact_pressure(n, d); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Int NTNBaseContact::getNodeIndex(UInt node) const { return this->slaves.find(node); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::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 << "NTNBaseContact [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + slaves : " << std::endl; this->slaves.printself(stream, indent + 2); stream << space << " + normals : " << std::endl; this->normals.printself(stream, indent + 2); stream << space << " + contact_pressure : " << std::endl; this->contact_pressure.printself(stream, indent + 2); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::syncArrays(SyncChoice sync_choice) { AKANTU_DEBUG_IN(); this->slaves.syncElements(sync_choice); this->normals.syncElements(sync_choice); this->is_in_contact.syncElements(sync_choice); this->contact_pressure.syncElements(sync_choice); this->lumped_boundary_slaves.syncElements(sync_choice); this->impedance.syncElements(sync_choice); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::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 dumper::NodalField, Array>( \ field, 0, 0, &nodal_filter)) if (field_id == "displacement") { ADD_FIELD(field_id, this->model.getDisplacement(), Real); } else if (field_id == "mass") { ADD_FIELD(field_id, this->model.getMass(), Real); } else if (field_id == "velocity") { ADD_FIELD(field_id, this->model.getVelocity(), Real); } else if (field_id == "acceleration") { ADD_FIELD(field_id, this->model.getAcceleration(), Real); } else if (field_id == "external_force") { ADD_FIELD(field_id, this->model.getForce(), Real); } else if (field_id == "internal_force") { ADD_FIELD(field_id, this->model.getInternalForce(), Real); } else if (field_id == "blocked_dofs") { ADD_FIELD(field_id, this->model.getBlockedDOFs(), bool); } else if (field_id == "increment") { ADD_FIELD(field_id, this->model.getIncrement(), Real); } else if (field_id == "normal") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->normals.getArray())); } else if (field_id == "contact_pressure") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->contact_pressure.getArray())); } else if (field_id == "is_in_contact") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->is_in_contact.getArray())); } else if (field_id == "lumped_boundary_slave") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->lumped_boundary_slaves.getArray())); } else if (field_id == "impedance") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->impedance.getArray())); } else { std::cerr << "Could not add field '" << field_id << "' to the dumper. Just ignored it." << std::endl; } #undef ADD_FIELD #endif AKANTU_DEBUG_OUT(); } } // namespace akantu 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 02b455c57..4e4232b23 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,587 +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 (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 (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 (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - Mesh::type_iterator it = mesh.firstType(dim - 1, *gt); - Mesh::type_iterator last = mesh.lastType(dim - 1, *gt); - - for (; it != last; ++it) { + 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, *it, - *gt); + boundary_fem.computeNormalsOnIntegrationPoints(cur_pos, quad_normals, type, + ghost_type); - UInt nb_quad_points = boundary_fem.getNbIntegrationPoints(*it, *gt); + 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(*it); - const Array & connectivity = mesh.getConnectivity(*it, *gt); + UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); + const Array & connectivity = mesh.getConnectivity(type, ghost_type); - Array::const_iterator elem_it = - (this->master_elements)(*it, *gt).begin(); - Array::const_iterator elem_it_end = - (this->master_elements)(*it, *gt).end(); // loop over contact nodes - for (; elem_it != elem_it_end; ++elem_it) { + for (auto & element : (this->master_elements)(type, ghost_type)) { for (UInt q = 0; q < nb_nodes_per_element; ++q) { - UInt node = connectivity(*elem_it, 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((*elem_it) * nb_quad_points + q, d); + quad_normals(element * nb_quad_points + q, d); } } } } - // ------------------------------------------------------------------------------------- - - /* - // get elements connected to each node - CSR node_to_element; - MeshUtils::buildNode2ElementsElementTypeMap(mesh, node_to_element, *it, - *gt); - - // add up normals to all master nodes - for (UInt n=0; nmasters(n); - CSR::iterator elem = node_to_element.begin(master); - // loop over all elements connected to this master node - for (; elem != node_to_element.end(master); ++elem) { - UInt e = *elem; - // loop over all quad points of this element - for (UInt q=0; qnormals(n,d) += quad_normals(e*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/packages/core.cmake b/packages/core.cmake index ffb02efa1..4ac17a8ca 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,522 +1,523 @@ #=============================================================================== # @file core.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Mon Jan 18 2016 # # @brief package description for core # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(core NOT_OPTIONAL DESCRIPTION "core package for Akantu" FEATURES_PUBLIC cxx_strong_enums cxx_defaulted_functions cxx_deleted_functions cxx_auto_type cxx_decltype_auto FEATURES_PRIVATE cxx_lambdas cxx_nullptr cxx_range_for cxx_delegating_constructors DEPENDS INTERFACE Boost) package_declare_sources(core common/aka_array.cc common/aka_array.hh common/aka_array_tmpl.hh common/aka_bbox.hh common/aka_blas_lapack.hh common/aka_circular_array.hh common/aka_circular_array_inline_impl.cc common/aka_common.cc common/aka_common.hh common/aka_common_inline_impl.cc common/aka_csr.hh common/aka_element_classes_info_inline_impl.cc + common/aka_enum_macros.hh common/aka_error.cc common/aka_error.hh common/aka_event_handler_manager.hh common/aka_extern.cc common/aka_factory.hh common/aka_fwd.hh common/aka_grid_dynamic.hh common/aka_math.cc common/aka_math.hh common/aka_math_tmpl.hh common/aka_memory.cc common/aka_memory.hh common/aka_memory_inline_impl.cc common/aka_named_argument.hh common/aka_random_generator.hh common/aka_safe_enum.hh common/aka_static_memory.cc common/aka_static_memory.hh common/aka_static_memory_inline_impl.cc common/aka_static_memory_tmpl.hh common/aka_typelist.hh common/aka_types.hh common/aka_visitor.hh common/aka_voigthelper.hh common/aka_voigthelper_tmpl.hh common/aka_voigthelper.cc common/aka_warning.hh common/aka_warning_restore.hh common/aka_iterators.hh common/aka_static_if.hh common/aka_compatibilty_with_cpp_standard.hh fe_engine/element_class.hh fe_engine/element_class_tmpl.hh fe_engine/element_classes/element_class_hexahedron_8_inline_impl.cc fe_engine/element_classes/element_class_hexahedron_20_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_6_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_15_inline_impl.cc fe_engine/element_classes/element_class_point_1_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_4_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc fe_engine/element_classes/element_class_segment_2_inline_impl.cc fe_engine/element_classes/element_class_segment_3_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.cc fe_engine/element_classes/element_class_triangle_3_inline_impl.cc fe_engine/element_classes/element_class_triangle_6_inline_impl.cc fe_engine/element_type_conversion.hh fe_engine/fe_engine.cc fe_engine/fe_engine.hh fe_engine/fe_engine_inline_impl.cc fe_engine/fe_engine_template.hh fe_engine/fe_engine_template_tmpl_field.hh fe_engine/fe_engine_template_tmpl.hh fe_engine/geometrical_element_property.hh fe_engine/geometrical_element_property.cc fe_engine/gauss_integration.cc fe_engine/gauss_integration_tmpl.hh fe_engine/integrator.hh fe_engine/integrator_gauss.hh fe_engine/integrator_gauss_inline_impl.cc fe_engine/interpolation_element_tmpl.hh fe_engine/integration_point.hh fe_engine/shape_functions.hh fe_engine/shape_functions.cc fe_engine/shape_functions_inline_impl.cc fe_engine/shape_lagrange_base.cc fe_engine/shape_lagrange_base.hh fe_engine/shape_lagrange_base_inline_impl.cc fe_engine/shape_lagrange.hh fe_engine/shape_lagrange_inline_impl.cc fe_engine/element.hh io/dumper/dumpable.hh io/dumper/dumpable.cc io/dumper/dumpable_dummy.hh io/dumper/dumpable_inline_impl.hh io/dumper/dumper_field.hh io/dumper/dumper_material_padders.hh io/dumper/dumper_filtered_connectivity.hh io/dumper/dumper_element_partition.hh io/mesh_io.cc io/mesh_io.hh io/mesh_io/mesh_io_abaqus.cc io/mesh_io/mesh_io_abaqus.hh io/mesh_io/mesh_io_diana.cc io/mesh_io/mesh_io_diana.hh io/mesh_io/mesh_io_msh.cc io/mesh_io/mesh_io_msh.hh #io/model_io.cc #io/model_io.hh io/parser/algebraic_parser.hh io/parser/input_file_parser.hh io/parser/parsable.cc io/parser/parsable.hh io/parser/parser.cc io/parser/parser_real.cc io/parser/parser_random.cc io/parser/parser_types.cc io/parser/parser_input_files.cc io/parser/parser.hh io/parser/parser_tmpl.hh io/parser/parser_grammar_tmpl.hh io/parser/cppargparse/cppargparse.hh io/parser/cppargparse/cppargparse.cc io/parser/cppargparse/cppargparse_tmpl.hh io/parser/parameter_registry.cc io/parser/parameter_registry.hh io/parser/parameter_registry_tmpl.hh mesh/element_group.cc mesh/element_group.hh mesh/element_group_inline_impl.cc mesh/element_type_map.cc mesh/element_type_map.hh mesh/element_type_map_tmpl.hh mesh/element_type_map_filter.hh mesh/group_manager.cc mesh/group_manager.hh mesh/group_manager_inline_impl.cc mesh/mesh.cc mesh/mesh.hh mesh/mesh_periodic.cc mesh/mesh_accessor.hh mesh/mesh_events.hh mesh/mesh_filter.hh mesh/mesh_global_data_updater.hh mesh/mesh_data.cc mesh/mesh_data.hh mesh/mesh_data_tmpl.hh mesh/mesh_inline_impl.cc mesh/node_group.cc mesh/node_group.hh mesh/node_group_inline_impl.cc mesh/mesh_iterators.hh mesh_utils/mesh_partition.cc mesh_utils/mesh_partition.hh mesh_utils/mesh_partition/mesh_partition_mesh_data.cc mesh_utils/mesh_partition/mesh_partition_mesh_data.hh mesh_utils/mesh_partition/mesh_partition_scotch.hh mesh_utils/mesh_utils_pbc.cc mesh_utils/mesh_utils.cc mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_distribution.cc mesh_utils/mesh_utils_distribution.hh mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_inline_impl.cc mesh_utils/global_ids_updater.hh mesh_utils/global_ids_updater.cc mesh_utils/global_ids_updater_inline_impl.cc model/boundary_condition.hh model/boundary_condition_functor.hh model/boundary_condition_functor_inline_impl.cc model/boundary_condition_tmpl.hh model/common/neighborhood_base.hh model/common/neighborhood_base.cc model/common/neighborhood_base_inline_impl.cc model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.hh model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion_inline_impl.cc model/common/non_local_toolbox/non_local_manager.hh model/common/non_local_toolbox/non_local_manager.cc model/common/non_local_toolbox/non_local_manager_inline_impl.cc model/common/non_local_toolbox/non_local_manager_callback.hh model/common/non_local_toolbox/non_local_neighborhood_base.hh model/common/non_local_toolbox/non_local_neighborhood_base.cc model/common/non_local_toolbox/non_local_neighborhood.hh model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh model/common/non_local_toolbox/non_local_neighborhood_inline_impl.cc model/common/non_local_toolbox/base_weight_function.hh model/common/non_local_toolbox/base_weight_function_inline_impl.cc model/dof_manager.cc model/dof_manager.hh model/dof_manager_default.cc model/dof_manager_default.hh model/dof_manager_default_inline_impl.cc model/dof_manager_inline_impl.cc model/model_solver.cc model/model_solver.hh model/non_linear_solver.cc model/non_linear_solver.hh model/non_linear_solver_default.hh model/non_linear_solver_lumped.cc model/non_linear_solver_lumped.hh model/solver_callback.hh model/solver_callback.cc model/time_step_solver.hh model/time_step_solvers/time_step_solver.cc model/time_step_solvers/time_step_solver_default.cc model/time_step_solvers/time_step_solver_default.hh model/time_step_solvers/time_step_solver_default_explicit.hh model/non_linear_solver_callback.hh model/time_step_solvers/time_step_solver_default_solver_callback.hh model/integration_scheme/generalized_trapezoidal.cc model/integration_scheme/generalized_trapezoidal.hh model/integration_scheme/integration_scheme.cc model/integration_scheme/integration_scheme.hh model/integration_scheme/integration_scheme_1st_order.cc model/integration_scheme/integration_scheme_1st_order.hh model/integration_scheme/integration_scheme_2nd_order.cc model/integration_scheme/integration_scheme_2nd_order.hh model/integration_scheme/newmark-beta.cc model/integration_scheme/newmark-beta.hh model/integration_scheme/pseudo_time.cc model/integration_scheme/pseudo_time.hh model/model.cc model/model.hh model/model_inline_impl.cc model/model_options.hh solver/sparse_solver.cc solver/sparse_solver.hh solver/sparse_solver_inline_impl.cc solver/sparse_matrix.cc solver/sparse_matrix.hh solver/sparse_matrix_inline_impl.cc solver/sparse_matrix_aij.cc solver/sparse_matrix_aij.hh solver/sparse_matrix_aij_inline_impl.cc solver/terms_to_assemble.hh synchronizer/communication_buffer_inline_impl.cc synchronizer/communication_descriptor.hh synchronizer/communication_descriptor_tmpl.hh synchronizer/communication_request.hh synchronizer/communication_tag.hh synchronizer/communications.hh synchronizer/communications_tmpl.hh synchronizer/communicator.cc synchronizer/communicator.hh synchronizer/communicator_dummy_inline_impl.cc synchronizer/communicator_event_handler.hh synchronizer/communicator_inline_impl.hh synchronizer/data_accessor.cc synchronizer/data_accessor.hh synchronizer/dof_synchronizer.cc synchronizer/dof_synchronizer.hh synchronizer/dof_synchronizer_inline_impl.cc synchronizer/element_info_per_processor.cc synchronizer/element_info_per_processor.hh synchronizer/element_info_per_processor_tmpl.hh synchronizer/element_synchronizer.cc synchronizer/element_synchronizer.hh synchronizer/facet_synchronizer.cc synchronizer/facet_synchronizer.hh synchronizer/facet_synchronizer_inline_impl.cc synchronizer/grid_synchronizer.cc synchronizer/grid_synchronizer.hh synchronizer/grid_synchronizer_tmpl.hh synchronizer/master_element_info_per_processor.cc synchronizer/node_info_per_processor.cc synchronizer/node_info_per_processor.hh synchronizer/node_synchronizer.cc synchronizer/node_synchronizer.hh synchronizer/periodic_node_synchronizer.cc synchronizer/periodic_node_synchronizer.hh synchronizer/slave_element_info_per_processor.cc synchronizer/synchronizer.cc synchronizer/synchronizer.hh synchronizer/synchronizer_impl.hh synchronizer/synchronizer_impl_tmpl.hh synchronizer/synchronizer_registry.cc synchronizer/synchronizer_registry.hh synchronizer/synchronizer_tmpl.hh synchronizer/communication_buffer.hh ) set(AKANTU_SPIRIT_SOURCES io/mesh_io/mesh_io_abaqus.cc io/parser/parser_real.cc io/parser/parser_random.cc io/parser/parser_types.cc io/parser/parser_input_files.cc PARENT_SCOPE ) package_declare_elements(core ELEMENT_TYPES _point_1 _segment_2 _segment_3 _triangle_3 _triangle_6 _quadrangle_4 _quadrangle_8 _tetrahedron_4 _tetrahedron_10 _pentahedron_6 _pentahedron_15 _hexahedron_8 _hexahedron_20 KIND regular GEOMETRICAL_TYPES _gt_point _gt_segment_2 _gt_segment_3 _gt_triangle_3 _gt_triangle_6 _gt_quadrangle_4 _gt_quadrangle_8 _gt_tetrahedron_4 _gt_tetrahedron_10 _gt_hexahedron_8 _gt_hexahedron_20 _gt_pentahedron_6 _gt_pentahedron_15 INTERPOLATION_TYPES _itp_lagrange_point_1 _itp_lagrange_segment_2 _itp_lagrange_segment_3 _itp_lagrange_triangle_3 _itp_lagrange_triangle_6 _itp_lagrange_quadrangle_4 _itp_serendip_quadrangle_8 _itp_lagrange_tetrahedron_4 _itp_lagrange_tetrahedron_10 _itp_lagrange_hexahedron_8 _itp_serendip_hexahedron_20 _itp_lagrange_pentahedron_6 _itp_lagrange_pentahedron_15 GEOMETRICAL_SHAPES _gst_point _gst_triangle _gst_square _gst_prism GAUSS_INTEGRATION_TYPES _git_point _git_segment _git_triangle _git_tetrahedron _git_pentahedron INTERPOLATION_KIND _itk_lagrangian FE_ENGINE_LISTS gradient_on_integration_points interpolate_on_integration_points interpolate compute_normals_on_integration_points inverse_map contains compute_shapes compute_shapes_derivatives get_shapes_derivatives lagrange_base ) package_declare_documentation_files(core manual.sty manual.cls manual.tex manual-macros.sty manual-titlepages.tex manual-authors.tex manual-changelog.tex manual-introduction.tex manual-gettingstarted.tex manual-io.tex manual-feengine.tex manual-elements.tex manual-appendix-elements.tex manual-appendix-packages.tex manual-backmatter.tex manual-bibliography.bib manual-bibliographystyle.bst figures/bc_and_ic_example.pdf figures/boundary.pdf figures/boundary.svg figures/dirichlet.pdf figures/dirichlet.svg # figures/doc_wheel.pdf # figures/doc_wheel.svg figures/hot-point-1.png figures/hot-point-2.png figures/insertion.pdf figures/interpolate.pdf figures/interpolate.svg figures/vectors.pdf figures/vectors.svg figures/elements/hexahedron_8.pdf figures/elements/hexahedron_8.svg figures/elements/quadrangle_4.pdf figures/elements/quadrangle_4.svg figures/elements/quadrangle_8.pdf figures/elements/quadrangle_8.svg figures/elements/segment_2.pdf figures/elements/segment_2.svg figures/elements/segment_3.pdf figures/elements/segment_3.svg figures/elements/tetrahedron_10.pdf figures/elements/tetrahedron_10.svg figures/elements/tetrahedron_4.pdf figures/elements/tetrahedron_4.svg figures/elements/triangle_3.pdf figures/elements/triangle_3.svg figures/elements/triangle_6.pdf figures/elements/triangle_6.svg figures/elements/xtemp.pdf ) package_declare_documentation(core "This package is the core engine of \\akantu. It depends on:" "\\begin{itemize}" "\\item A C++ compiler (\\href{http://gcc.gnu.org/}{GCC} >= 4, or \\href{https://software.intel.com/en-us/intel-compilers}{Intel})." "\\item The cross-platform, open-source \\href{http://www.cmake.org/}{CMake} build system." "\\item The \\href{http://www.boost.org/}{Boost} C++ portable libraries." "\\item The \\href{http://www.zlib.net/}{zlib} compression library." "\\end{itemize}" "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "\\begin{command}" " > sudo apt-get install cmake libboost-dev zlib1g-dev g++" "\\end{command}" "" "Under Mac OS X the installation requires the following steps:" "\\begin{itemize}" "\\item Install Xcode" "\\item Install the command line tools." "\\item Install the MacPorts project which allows to automatically" "download and install opensource packages." "\\end{itemize}" "Then the following commands should be typed in a terminal:" "\\begin{command}" " > sudo port install cmake gcc48 boost" "\\end{command}" ) find_program(READLINK_COMMAND readlink) find_program(ADDR2LINE_COMMAND addr2line) find_program(PATCH_COMMAND patch) mark_as_advanced(READLINK_COMMAND) mark_as_advanced(ADDR2LINE_COMMAND) package_declare_extra_files_to_package(core SOURCES common/aka_element_classes_info.hh.in common/aka_config.hh.in ) if(CMAKE_CXX_COMPILER_ID STREQUAL "Clang" AND (NOT CMAKE_CXX_COMPILER_VERSION VERSION_LESS 3.9)) package_set_compile_flags(core CXX "-Wno-undefined-var-template") endif() if(DEFINED AKANTU_CXX11_FLAGS) package_declare(core_cxx11 NOT_OPTIONAL DESCRIPTION "C++ 11 additions for Akantu core" COMPILE_FLAGS CXX "${AKANTU_CXX11_FLAGS}") if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "4.6") set(AKANTU_CORE_CXX11 OFF CACHE BOOL "C++ 11 additions for Akantu core - not supported by the selected compiler" FORCE) endif() endif() package_declare_documentation(core_cxx11 "This option activates some features of the C++11 standard. This is usable with GCC>=4.7 or Intel>=13.") else() if(CMAKE_VERSION VERSION_LESS 3.1) message(FATAL_ERROR "Since version 3.0 Akantu requires at least c++11 capable compiler") endif() endif() diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 63d808f9f..1ece02405 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,564 +1,493 @@ /** * @file aka_common.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Feb 12 2018 * * @brief common type descriptions for akantu * * @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 * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ #include "aka_compatibilty_with_cpp_standard.hh" /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU_DUMPER__ namespace dumper { #define __END_AKANTU_DUMPER__ } /* -------------------------------------------------------------------------- */ #if defined(WIN32) #define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ using ID = std::string; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = Real(0.); #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ using MemoryID = UInt; // using Surface = std::string; // using SurfacePair= std::pair; // using SurfacePairList = std::list; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; -#define AKANTU_PP_ENUM(s, data, i, elem) \ - BOOST_PP_TUPLE_REM() \ - elem BOOST_PP_COMMA_IF(BOOST_PP_NOT_EQUAL(i, BOOST_PP_DEC(data))) -} // namespace akantu - -#if (defined(__GNUC__) || defined(__GNUG__)) -#define AKA_GCC_VERSION \ - (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) -#if AKA_GCC_VERSION < 60000 -#define AKANTU_ENUM_HASH(type_name) \ - namespace std { \ - template <> struct hash<::akantu::type_name> { \ - using argument_type = ::akantu::type_name; \ - size_t operator()(const argument_type & e) const noexcept { \ - auto ue = underlying_type_t(e); \ - return uh(ue); \ - } \ - \ - private: \ - const hash> uh{}; \ - }; \ - } -#else -#define AKANTU_ENUM_HASH(type_name) -#endif // AKA_GCC_VERSION -#endif // GNU - +#include "aka_enum_macros.hh" +/* -------------------------------------------------------------------------- */ #include "aka_element_classes_info.hh" namespace akantu { - -#define AKANTU_PP_CAT(s, data, elem) BOOST_PP_CAT(data, elem) - -#define AKANTU_PP_TYPE_TO_STR(s, data, elem) \ - ({BOOST_PP_CAT(data::_, elem), BOOST_PP_STRINGIZE(elem)}) - -#define AKANTU_PP_STR_TO_TYPE(s, data, elem) \ - ({BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data::_, elem)}) - -#define AKANTU_ENUM_DECLARE(type_name, list) \ - enum class type_name { \ - BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_CAT, _, list)) \ - }; - -#define AKANTU_ENUM_OUTPUT_STREAM(type_name, list) \ - } \ - AKANTU_ENUM_HASH(type_name) \ - namespace aka { \ - inline std::string to_string(const ::akantu::type_name & type) { \ - static std::unordered_map<::akantu::type_name, std::string> convert{ \ - BOOST_PP_SEQ_FOR_EACH_I( \ - AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ - BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_TYPE_TO_STR, \ - ::akantu::type_name, list))}; \ - return convert.at(type); \ - } \ - } \ - namespace akantu { \ - inline std::ostream & operator<<(std::ostream & stream, \ - const type_name & type) { \ - stream << aka::to_string(type); \ - return stream; \ - } - -#define AKANTU_ENUM_INPUT_STREAM(type_name, list) \ - inline std::istream & operator>>(std::istream & stream, type_name & type) { \ - std::string str; \ - stream >> str; \ - static std::unordered_map convert{ \ - BOOST_PP_SEQ_FOR_EACH_I( \ - AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(list), \ - BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE, type_name, list))}; \ - type = convert.at(str); \ - return stream; \ - } - /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ /// small help to use names for directions enum SpatialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has /// structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum MeshEventHandlerPriority defines relative order of execution of events enum EventHandlerPriority { _ehp_highest = 0, _ehp_mesh = 5, _ehp_fe_engine = 9, _ehp_synchronizer = 10, _ehp_dof_manager = 20, _ehp_model = 94, _ehp_non_local_manager = 100, _ehp_lowest = 100 }; #ifndef SWIG // clang-format off #define AKANTU_MODEL_TYPES \ (model) \ (solid_mechanics_model) \ (solid_mechanics_model_cohesive) \ (heat_transfer_model) \ (structural_mechanics_model) \ (embedded_model) // clang-format on /// enum ModelType defines which type of physics is solved -AKANTU_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) -AKANTU_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) -AKANTU_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) +AKANTU_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES) +AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) +AKANTU_CLASS_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES) #else enum class ModelType { _model, _solid_mechanics_model, _solid_mechanics_model_cohesive, _heat_transfer_model, _structural_mechanics_model, _embedded_model }; #endif /// enum AnalysisMethod type of solving method used to solve the equation of /// motion enum AnalysisMethod { _static = 0, _implicit_dynamic = 1, _explicit_lumped_mass = 2, _explicit_lumped_capacity = 2, _explicit_consistent_mass = 3 }; /// enum DOFSupportType defines which kind of dof that can exists enum DOFSupportType { _dst_nodal, _dst_generic }; /// Type of non linear resolution available in akantu enum NonLinearSolverType { _nls_linear, ///< No non linear convergence loop _nls_newton_raphson, ///< Regular Newton-Raphson _nls_newton_raphson_modified, ///< Newton-Raphson with initial tangent _nls_lumped, ///< Case of lumped mass or equivalent matrix _nls_auto ///< This will take a default value that make sense in case of /// model::getNewSolver }; /// Type of time stepping solver enum TimeStepSolverType { _tsst_static, ///< Static solution _tsst_dynamic, ///< Dynamic solver _tsst_dynamic_lumped, ///< Dynamic solver with lumped mass _tsst_not_defined, ///< For not defined cases }; /// Type of integration scheme enum IntegrationSchemeType { _ist_pseudo_time, ///< Pseudo Time _ist_forward_euler, ///< GeneralizedTrapezoidal(0) _ist_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2) _ist_backward_euler, ///< GeneralizedTrapezoidal(1) _ist_central_difference, ///< NewmarkBeta(0, 1/2) _ist_fox_goodwin, ///< NewmarkBeta(1/6, 1/2) _ist_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2) _ist_linear_acceleration, ///< NewmarkBeta(1/3, 1/2) _ist_newmark_beta, ///< generic NewmarkBeta with user defined /// alpha and beta _ist_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user /// defined alpha }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_solution, ///< Use solution to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// @enum SparseMatrixType type of sparse matrix used enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- Generic tags --- _gst_whatever, _gst_update, _gst_ask_nodes, _gst_size, //--- SolidMechanicsModel tags --- _gst_smm_mass, ///< synchronization of the SolidMechanicsModel.mass _gst_smm_for_gradu, ///< synchronization of the /// SolidMechanicsModel.displacement _gst_smm_boundary, ///< synchronization of the boundary, forces, velocities /// and displacement _gst_smm_uv, ///< synchronization of the nodal velocities and displacement _gst_smm_res, ///< synchronization of the nodal residual _gst_smm_init_mat, ///< synchronization of the data to initialize materials _gst_smm_stress, ///< synchronization of the stresses to compute the internal /// forces _gst_smmc_facets, ///< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, ///< synchronization of facet global connectivity _gst_smmc_facets_stress, ///< synchronization of facets' stress to setup facet /// synch _gst_smmc_damage, ///< synchronization of damage // --- GlobalIdsUpdater tags --- _gst_giu_global_conn, ///< synchronization of global connectivities // --- CohesiveElementInserter tags --- _gst_ce_groups, ///< synchronization of cohesive element insertion depending /// on facet groups // --- GroupManager tags --- _gst_gm_clusters, ///< synchronization of clusters // --- HeatTransfer tags --- _gst_htm_temperature, ///< synchronization of the nodal temperature _gst_htm_gradient_temperature, ///< synchronization of the element gradient /// temperature // --- LevelSet tags --- _gst_htm_phi, ///< synchronization of the nodal level set value phi _gst_htm_gradient_phi, ///< synchronization of the element gradient phi //--- Material non local --- _gst_mnl_for_average, ///< synchronization of data to average in non local /// material _gst_mnl_weight, ///< synchronization of data for the weight computations // --- NeighborhoodSynchronization tags --- _gst_nh_criterion, // --- General tags --- _gst_test, ///< Test tag _gst_user_1, ///< tag for user simulations _gst_user_2, ///< tag for user simulations _gst_material_id, ///< synchronization of the material ids _gst_for_dump, ///< everything that needs to be synch before dump // --- Contact & Friction --- _gst_cf_nodal, ///< synchronization of disp, velo, and current position _gst_cf_incr, ///< synchronization of increment // --- Solver tags --- _gst_solver_solution ///< synchronization of the solution obained with the /// PETSc solver }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator<<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost = 0, _ghost = 1, _casper // not used but a real cute ghost }; /// Define the flag that can be set to a node enum class NodeFlag : std::uint8_t { _normal = 0x00, _distributed = 0x01, _master = 0x03, _slave = 0x05, _pure_ghost = 0x09, _shared_mask = 0x0F, _periodic = 0x10, _periodic_master = 0x30, _periodic_slave = 0x50, _periodic_mask = 0xF0, _local_master_mask = 0xCC, // ~(_master & _periodic_mask) }; inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) & under(b)); } inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) { using under = std::underlying_type_t; return NodeFlag(under(a) | under(b)); } inline NodeFlag & operator|=(NodeFlag & a, const NodeFlag & b) { a = a | b; return a; } inline NodeFlag & operator&=(NodeFlag & a, const NodeFlag & b) { a = a & b; return a; } inline NodeFlag operator~(const NodeFlag & a) { using under = std::underlying_type_t; return NodeFlag(~under(a)); } inline std::ostream & operator<<(std::ostream & stream, const NodeFlag & flag) { using under = std::underlying_type_t; stream << under(flag); return stream; } } // namespace akantu #ifndef SWIG AKANTU_ENUM_HASH(GhostType) #endif namespace akantu { /* -------------------------------------------------------------------------- */ struct GhostType_def { using type = GhostType; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; using ghost_type_t = safe_enum; extern ghost_type_t ghost_types; /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type); /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ /* Type traits */ /* -------------------------------------------------------------------------- */ struct TensorTrait {}; /* -------------------------------------------------------------------------- */ template using is_tensor = std::is_base_of; /* -------------------------------------------------------------------------- */ template using is_scalar = std::is_arithmetic; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name(type variable) { this->variable = variable; } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name() const { return variable; } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name() { return variable; } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \ inline con Array & get##name( \ const support & el_type, const GhostType & ghost_type = _not_ghost) \ con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, ) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, ) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char **& argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char **& argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize(); /* -------------------------------------------------------------------------- */ /// Read an new input file void readInputFile(const std::string & input_file); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ // #if defined(__INTEL_COMPILER) // /// remark #981: operands are evaluated in unspecified order // #pragma warning(disable : 981) // /// remark #383: value copied to temporary, reference to temporary used // #pragma warning(disable : 383) // #endif // defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); inline std::string trim(const std::string & to_trim, char c); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ } // namespace akantu #include "aka_fwd.hh" namespace akantu { /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); } // namespace akantu #include "aka_common_inline_impl.cc" /* -------------------------------------------------------------------------- */ #if AKANTU_INTEGER_SIZE == 4 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9 #elif AKANTU_INTEGER_SIZE == 8 #define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL #endif namespace std { /** * Hashing function for pairs based on hash_combine from boost The magic number * is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f] * @f[\frac{2^32}{\phi} = 0x9e3779b9@f] * http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine * http://burtleburtle.net/bob/hash/doobs.html */ template struct hash> { hash() = default; size_t operator()(const std::pair & p) const { size_t seed = ah(p.first); return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) + (seed >> 2); } private: const hash ah{}; const hash bh{}; }; } // namespace std #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/common/aka_element_classes_info_inline_impl.cc b/src/common/aka_element_classes_info_inline_impl.cc index 871c7d2b6..bcf9cb8ea 100644 --- a/src/common/aka_element_classes_info_inline_impl.cc +++ b/src/common/aka_element_classes_info_inline_impl.cc @@ -1,114 +1,53 @@ /** * @file aka_element_classes_info_inline_impl.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Thu Jun 18 2015 * @date last modification: Wed Jan 10 2018 * * @brief Implementation of the streaming fonction for the element classes * enums * * @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 /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_ELEMENT_CLASSES_INFO_INLINE_IMPL_CC__ #define __AKANTU_AKA_ELEMENT_CLASSES_INFO_INLINE_IMPL_CC__ -AKANTU_ENUM_HASH(ElementType) - -#define AKANTU_PP_ELEMTYPE_TO_STR(s, data, elem) \ - ({akantu::elem, BOOST_PP_STRINGIZE(elem)}) - -#define AKANTU_PP_STR_TO_ELEMTYPE(s, data, elem) \ - ({BOOST_PP_STRINGIZE(elem), akantu::elem}) - -namespace aka { -inline std::string to_string(const akantu::ElementType & type) { - static std::unordered_map convert{ - BOOST_PP_SEQ_FOR_EACH_I( - AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(AKANTU_ALL_ELEMENT_TYPE), - BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_ELEMTYPE_TO_STR, _, - AKANTU_ALL_ELEMENT_TYPE)), - {akantu::_not_defined, "_not_defined"}, - {akantu::_max_element_type, "_max_element_type"}}; - return convert.at(type); -} -} - namespace akantu { -#define STRINGIFY(type) stream << BOOST_PP_STRINGIZE(type) +AKANTU_ENUM_OUTPUT_STREAM(ElementType, AKANTU_ALL_ELEMENT_TYPE(_not_defined)(_max_element_type)) +AKANTU_ENUM_INPUT_STREAM(ElementType, AKANTU_ALL_ELEMENT_TYPE) -/* -------------------------------------------------------------------------- */ -//! standard output stream operator for ElementType -inline std::ostream & operator<<(std::ostream & stream, - const ElementType & type) { - stream << aka::to_string(type); - return stream; -} - -/* -------------------------------------------------------------------------- */ - -//! standard input stream operator for ElementType -inline std::istream & operator>>(std::istream & stream, ElementType & type) { - std::string str; - stream >> str; - static std::unordered_map convert{ - BOOST_PP_SEQ_FOR_EACH_I( - AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(AKANTU_ALL_ELEMENT_TYPE), - BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_ELEMTYPE, _, - AKANTU_ALL_ELEMENT_TYPE))}; - type = convert.at(str); - return stream; -} - -/* -------------------------------------------------------------------------- */ -//! standard output stream operator for ElementType -inline std::ostream & operator<<(std::ostream & stream, ElementKind kind) { - AKANTU_BOOST_ALL_KIND_SWITCH(STRINGIFY); - - return stream; -} - -/* -------------------------------------------------------------------------- */ -/// standard output stream operator for InterpolationType -inline std::ostream & operator<<(std::ostream & stream, - InterpolationType type) { - switch (type) { - BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, STRINGIFY, - AKANTU_INTERPOLATION_TYPES) - case _itp_not_defined: - stream << "_itp_not_defined"; - break; - } - return stream; -} +AKANTU_ENUM_OUTPUT_STREAM(InterpolationType, AKANTU_INTERPOLATION_TYPES) +AKANTU_ENUM_INPUT_STREAM(InterpolationType, AKANTU_INTERPOLATION_TYPES) -#undef STRINGIFY +AKANTU_ENUM_OUTPUT_STREAM(ElementKind, AKANTU_ELEMENT_KIND) +AKANTU_ENUM_INPUT_STREAM(ElementKind, AKANTU_ELEMENT_KIND) -} // akantu +} // namespace akantu #endif /* __AKANTU_AKA_ELEMENT_CLASSES_INFO_INLINE_IMPL_CC__ */ diff --git a/src/fe_engine/shape_functions.cc b/src/fe_engine/shape_functions.cc index b45c222a9..ae1eebadc 100644 --- a/src/fe_engine/shape_functions.cc +++ b/src/fe_engine/shape_functions.cc @@ -1,245 +1,233 @@ /** * @file shape_functions.cc * * @author Nicolas Richart * * @date creation: Wed Aug 09 2017 * @date last modification: Wed Oct 11 2017 * * @brief implementation of th shape functions interface * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "shape_functions.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ShapeFunctions::ShapeFunctions(const Mesh & mesh, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), shapes("shapes_generic", id, memory_id), shapes_derivatives("shapes_derivatives_generic", id, memory_id), mesh(mesh) {} /* -------------------------------------------------------------------------- */ template inline void ShapeFunctions::initElementalFieldInterpolationFromIntegrationPoints( const Array & interpolation_points_coordinates, ElementTypeMapArray & interpolation_points_coordinates_matrices, ElementTypeMapArray & quad_points_coordinates_inv_matrices, const Array & quadrature_points_coordinates, const GhostType & ghost_type, const Array & element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); UInt nb_element = this->mesh.getNbElement(type, ghost_type); UInt nb_element_filter; if (element_filter == empty_filter) nb_element_filter = nb_element; else nb_element_filter = element_filter.size(); UInt nb_quad_per_element = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_interpolation_points_per_elem = interpolation_points_coordinates.size() / nb_element; AKANTU_DEBUG_ASSERT(interpolation_points_coordinates.size() % nb_element == 0, "Number of interpolation points should be a multiple of " "total number of elements"); if (!quad_points_coordinates_inv_matrices.exists(type, ghost_type)) quad_points_coordinates_inv_matrices.alloc( nb_element_filter, nb_quad_per_element * nb_quad_per_element, type, ghost_type); else quad_points_coordinates_inv_matrices(type, ghost_type) .resize(nb_element_filter); if (!interpolation_points_coordinates_matrices.exists(type, ghost_type)) interpolation_points_coordinates_matrices.alloc( nb_element_filter, nb_interpolation_points_per_elem * nb_quad_per_element, type, ghost_type); else interpolation_points_coordinates_matrices(type, ghost_type) .resize(nb_element_filter); Array & quad_inv_mat = quad_points_coordinates_inv_matrices(type, ghost_type); Array & interp_points_mat = interpolation_points_coordinates_matrices(type, ghost_type); Matrix quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); Array::const_matrix_iterator quad_coords_it = quadrature_points_coordinates.begin_reinterpret( spatial_dimension, nb_quad_per_element, nb_element_filter); Array::const_matrix_iterator points_coords_begin = interpolation_points_coordinates.begin_reinterpret( spatial_dimension, nb_interpolation_points_per_elem, nb_element); Array::matrix_iterator inv_quad_coord_it = quad_inv_mat.begin(nb_quad_per_element, nb_quad_per_element); Array::matrix_iterator int_points_mat_it = interp_points_mat.begin( nb_interpolation_points_per_elem, nb_quad_per_element); /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element_filter; ++el, ++inv_quad_coord_it, ++int_points_mat_it, ++quad_coords_it) { /// matrix containing the quadrature points coordinates const Matrix & quad_coords = *quad_coords_it; /// matrix to store the matrix inversion result Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /// insert the quad coordinates in a matrix compatible with the /// interpolation buildElementalFieldInterpolationMatrix(quad_coords, quad_coord_matrix); /// invert the interpolation matrix inv_quad_coord_matrix.inverse(quad_coord_matrix); /// matrix containing the interpolation points coordinates const Matrix & points_coords = points_coords_begin[element_filter(el)]; /// matrix to store the interpolation points coordinates /// compatible with these functions Matrix & inv_points_coord_matrix = *int_points_mat_it; /// insert the quad coordinates in a matrix compatible with the /// interpolation buildElementalFieldInterpolationMatrix(points_coords, inv_points_coord_matrix); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ShapeFunctions::initElementalFieldInterpolationFromIntegrationPoints( const ElementTypeMapArray & interpolation_points_coordinates, ElementTypeMapArray & interpolation_points_coordinates_matrices, ElementTypeMapArray & quad_points_coordinates_inv_matrices, const ElementTypeMapArray & quadrature_points_coordinates, const ElementTypeMapArray * element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); for (auto ghost_type : ghost_types) { - Mesh::type_iterator it, last; + auto types_iterable = mesh.elementTypes(spatial_dimension, ghost_type); if (element_filter) { - it = element_filter->firstType(spatial_dimension, ghost_type); - last = element_filter->lastType(spatial_dimension, ghost_type); - } else { - it = mesh.firstType(spatial_dimension, ghost_type); - last = mesh.lastType(spatial_dimension, ghost_type); + types_iterable = element_filter->elementTypes(spatial_dimension, ghost_type); } - for (; it != last; ++it) { - ElementType type = *it; + for (auto type : types_iterable) { UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; const Array * elem_filter; if (element_filter) elem_filter = &((*element_filter)(type, ghost_type)); else elem_filter = &(empty_filter); #define AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS(type) \ this->initElementalFieldInterpolationFromIntegrationPoints( \ interpolation_points_coordinates(type, ghost_type), \ interpolation_points_coordinates_matrices, \ quad_points_coordinates_inv_matrices, \ quadrature_points_coordinates(type, ghost_type), ghost_type, \ *elem_filter) AKANTU_BOOST_REGULAR_ELEMENT_SWITCH( AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS); #undef AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ShapeFunctions::interpolateElementalFieldFromIntegrationPoints( const ElementTypeMapArray & field, const ElementTypeMapArray & interpolation_points_coordinates_matrices, const ElementTypeMapArray & quad_points_coordinates_inv_matrices, ElementTypeMapArray & result, const GhostType & ghost_type, const ElementTypeMapArray * element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); - Mesh::type_iterator it, last; - + auto types_iterable = mesh.elementTypes(spatial_dimension, ghost_type); if (element_filter) { - it = element_filter->firstType(spatial_dimension, ghost_type); - last = element_filter->lastType(spatial_dimension, ghost_type); - } else { - it = mesh.firstType(spatial_dimension, ghost_type); - last = mesh.lastType(spatial_dimension, ghost_type); + types_iterable = element_filter->elementTypes(spatial_dimension, ghost_type); } - for (; it != last; ++it) { - - ElementType type = *it; + for (auto type : types_iterable) { UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; const Array * elem_filter; if (element_filter) elem_filter = &((*element_filter)(type, ghost_type)); else elem_filter = &(empty_filter); #define AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS(type) \ interpolateElementalFieldFromIntegrationPoints( \ field(type, ghost_type), \ interpolation_points_coordinates_matrices(type, ghost_type), \ quad_points_coordinates_inv_matrices(type, ghost_type), result, \ ghost_type, *elem_filter) AKANTU_BOOST_REGULAR_ELEMENT_SWITCH( AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS); #undef AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/io/dumper/dumper_compute.hh b/src/io/dumper/dumper_compute.hh index 0989f62b2..f16c1c9e5 100644 --- a/src/io/dumper/dumper_compute.hh +++ b/src/io/dumper/dumper_compute.hh @@ -1,274 +1,268 @@ /** * @file dumper_compute.hh * * @author Guillaume Anciaux * * @date creation: Tue Sep 02 2014 * @date last modification: Sun Dec 03 2017 * * @brief Field that map a function to another field * * @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_DUMPER_COMPUTE_HH__ #define __AKANTU_DUMPER_COMPUTE_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "dumper_field.hh" #include "dumper_iohelper.hh" #include "dumper_type_traits.hh" #include /* -------------------------------------------------------------------------- */ namespace akantu { __BEGIN_AKANTU_DUMPER__ class ComputeFunctorInterface { public: virtual ~ComputeFunctorInterface() = default; virtual UInt getDim() = 0; virtual UInt getNbComponent(UInt old_nb_comp) = 0; }; /* -------------------------------------------------------------------------- */ template class ComputeFunctorOutput : public ComputeFunctorInterface { public: ComputeFunctorOutput() = default; ~ComputeFunctorOutput() override = default; }; /* -------------------------------------------------------------------------- */ template class ComputeFunctor : public ComputeFunctorOutput { public: ComputeFunctor() = default; ~ComputeFunctor() override = default; virtual return_type func(const input_type & d, Element global_index) = 0; }; /* -------------------------------------------------------------------------- */ template class FieldCompute : public Field { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ public: using sub_iterator = typename SubFieldCompute::iterator; using sub_types = typename SubFieldCompute::types; using sub_return_type = typename sub_types::return_type; using return_type = _return_type; using data_type = typename sub_types::data_type; using types = TypeTraits>; class iterator { public: iterator(const sub_iterator & it, ComputeFunctor & func) : it(it), func(func) {} bool operator!=(const iterator & it) const { return it.it != this->it; } iterator operator++() { ++this->it; return *this; } UInt currentGlobalIndex() { return this->it.currentGlobalIndex(); } return_type operator*() { return func.func(*it, it.getCurrentElement()); } Element getCurrentElement() { return this->it.getCurrentElement(); } UInt element_type() { return this->it.element_type(); } protected: sub_iterator it; ComputeFunctor & func; }; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FieldCompute(SubFieldCompute & cont, ComputeFunctorInterface & func) : sub_field(cont), func(dynamic_cast &>( func)) { this->checkHomogeneity(); }; ~FieldCompute() override { delete &(this->sub_field); delete &(this->func); } void registerToDumper(const std::string & id, iohelper::Dumper & dumper) override { dumper.addElemDataField(id, *this); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: iterator begin() { return iterator(sub_field.begin(), func); } iterator end() { return iterator(sub_field.end(), func); } UInt getDim() { return func.getDim(); } UInt size() { throw; // return Functor::size(); return 0; } void checkHomogeneity() override { this->homogeneous = true; }; iohelper::DataType getDataType() { return iohelper::getDataType(); } /// get the number of components of the hosted field ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) override { ElementTypeMap nb_components; - const ElementTypeMap & old_nb_components = + const auto & old_nb_components = this->sub_field.getNbComponents(dim, ghost_type, kind); - ElementTypeMap::type_iterator tit = - old_nb_components.firstType(dim, ghost_type, kind); - ElementTypeMap::type_iterator end = - old_nb_components.lastType(dim, ghost_type, kind); - - while (tit != end) { - UInt nb_comp = old_nb_components(*tit, ghost_type); - nb_components(*tit, ghost_type) = func.getNbComponent(nb_comp); - ++tit; + for(auto type : old_nb_components.elementTypes(dim, ghost_type, kind)) { + UInt nb_comp = old_nb_components(type, ghost_type); + nb_components(type, ghost_type) = func.getNbComponent(nb_comp); } return nb_components; }; /// for connection to a FieldCompute inline Field * connect(FieldComputeProxy & proxy) override; /// for connection to a FieldCompute ComputeFunctorInterface * connect(HomogenizerProxy & proxy) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: SubFieldCompute & sub_field; ComputeFunctor & func; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class FieldComputeProxy { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FieldComputeProxy(ComputeFunctorInterface & func) : func(func){}; inline static Field * createFieldCompute(Field * field, ComputeFunctorInterface & func) { /// that looks fishy an object passed as a ref and destroyed at their end of /// the function FieldComputeProxy compute_proxy(func); return field->connect(compute_proxy); } template Field * connectToField(T * ptr) { if (dynamic_cast> *>(&func)) { return this->connectToFunctor>(ptr); } else if (dynamic_cast> *>(&func)) { return this->connectToFunctor>(ptr); } else if (dynamic_cast> *>(&func)) { return this->connectToFunctor>(ptr); } else if (dynamic_cast> *>(&func)) { return this->connectToFunctor>(ptr); } else throw; } template Field * connectToFunctor(T * ptr) { auto * functor_ptr = new FieldCompute(*ptr, func); return functor_ptr; } template Field * connectToFunctor(__attribute__((unused)) FieldCompute, return_type2> * ptr) { throw; // return new FieldCompute(*ptr,func); return nullptr; } template Field * connectToFunctor( __attribute__((unused)) FieldCompute< FieldCompute, return_type2>, return_type3>, return_type4> * ptr) { throw; // return new FieldCompute(*ptr,func); return nullptr; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: ComputeFunctorInterface & func; }; /* -------------------------------------------------------------------------- */ /// for connection to a FieldCompute template inline Field * FieldCompute::connect(FieldComputeProxy & proxy) { return proxy.connectToField(this); } /* -------------------------------------------------------------------------- */ __END_AKANTU_DUMPER__ } // akantu #endif /* __AKANTU_DUMPER_COMPUTE_HH__ */ diff --git a/src/io/dumper/dumper_generic_elemental_field.hh b/src/io/dumper/dumper_generic_elemental_field.hh index f616aa8b7..a66095318 100644 --- a/src/io/dumper/dumper_generic_elemental_field.hh +++ b/src/io/dumper/dumper_generic_elemental_field.hh @@ -1,225 +1,215 @@ /** * @file dumper_generic_elemental_field.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Sep 02 2014 * @date last modification: Wed Nov 08 2017 * * @brief Generic interface for elemental fields * * @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_DUMPER_GENERIC_ELEMENTAL_FIELD_HH__ #define __AKANTU_DUMPER_GENERIC_ELEMENTAL_FIELD_HH__ /* -------------------------------------------------------------------------- */ #include "dumper_element_iterator.hh" #include "dumper_field.hh" #include "dumper_homogenizing_field.hh" #include "element_type_map_filter.hh" /* -------------------------------------------------------------------------- */ namespace akantu { __BEGIN_AKANTU_DUMPER__ /* -------------------------------------------------------------------------- */ template class iterator_type> class GenericElementalField : public Field { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ public: // check dumper_type_traits.hh for additional information over these types using types = _types; using data_type = typename types::data_type; using it_type = typename types::it_type; using field_type = typename types::field_type; using array_type = typename types::array_type; using array_iterator = typename types::array_iterator; using field_type_iterator = typename field_type::type_iterator; using iterator = iterator_type; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GenericElementalField(const field_type & field, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) : field(field), spatial_dimension(spatial_dimension), ghost_type(ghost_type), element_kind(element_kind) { this->checkHomogeneity(); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get the number of components of the hosted field ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) override { return this->field.getNbComponents(dim, ghost_type, kind); }; /// return the size of the contained data: i.e. the number of elements ? virtual UInt size() { checkHomogeneity(); return this->nb_total_element; } /// return the iohelper datatype to be dumped iohelper::DataType getDataType() { return iohelper::getDataType(); } protected: /// return the number of entries per element UInt getNbDataPerElem(const ElementType & type, const GhostType & ghost_type = _not_ghost) const { if (!nb_data_per_elem.exists(type, ghost_type)) return field(type, ghost_type).getNbComponent(); return nb_data_per_elem(type, this->ghost_type); } /// check if the same quantity of data for all element types void checkHomogeneity() override; public: void registerToDumper(const std::string & id, iohelper::Dumper & dumper) override { dumper.addElemDataField(id, *this); }; /// for connection to a FieldCompute inline Field * connect(FieldComputeProxy & proxy) override { return proxy.connectToField(this); } /// for connection to a Homogenizer inline ComputeFunctorInterface * connect(HomogenizerProxy & proxy) override { return proxy.connectToField(this); }; virtual iterator begin() { - field_type_iterator tit; - field_type_iterator end; - /// type iterators on the elemental field - tit = this->field.firstType(this->spatial_dimension, this->ghost_type, - this->element_kind); - end = this->field.lastType(this->spatial_dimension, this->ghost_type, - this->element_kind); + auto types = this->field.elementTypes(this->spatial_dimension, this->ghost_type, + this->element_kind); + auto tit = types.begin(); + auto end = types.end(); /// skip all types without data - ElementType type = *tit; for (; tit != end && this->field(*tit, this->ghost_type).size() == 0; ++tit) { } - type = *tit; + auto type = *tit; if (tit == end) return this->end(); /// getting information for the field of the given type - const array_type & vect = this->field(type, this->ghost_type); + const auto & vect = this->field(type, this->ghost_type); UInt nb_data_per_elem = this->getNbDataPerElem(type); - UInt nb_component = vect.getNbComponent(); - UInt size = (vect.size() * nb_component) / nb_data_per_elem; /// define element-wise iterator - array_iterator it = vect.begin_reinterpret(nb_data_per_elem, size); - array_iterator it_end = vect.end_reinterpret(nb_data_per_elem, size); + auto view = make_view(vect, nb_data_per_elem); + auto it = view.begin(); + auto it_end = view.end(); /// define data iterator iterator rit = iterator(this->field, tit, end, it, it_end, this->ghost_type); rit.setNbDataPerElem(this->nb_data_per_elem); return rit; } virtual iterator end() { - field_type_iterator tit; - field_type_iterator end; - - tit = this->field.firstType(this->spatial_dimension, this->ghost_type, - this->element_kind); - end = this->field.lastType(this->spatial_dimension, this->ghost_type, - this->element_kind); + auto types = this->field.elementTypes(this->spatial_dimension, this->ghost_type, + this->element_kind); + auto tit = types.begin(); + auto end = types.end(); - ElementType type = *tit; + auto type = *tit; for (; tit != end; ++tit) type = *tit; const array_type & vect = this->field(type, this->ghost_type); UInt nb_data = this->getNbDataPerElem(type); - UInt nb_component = vect.getNbComponent(); - UInt size = (vect.size() * nb_component) / nb_data; - array_iterator it = vect.end_reinterpret(nb_data, size); - - iterator rit = iterator(this->field, end, end, it, it, this->ghost_type); + auto it = make_view(vect, nb_data).end(); + auto rit = iterator(this->field, end, end, it, it, this->ghost_type); rit.setNbDataPerElem(this->nb_data_per_elem); + return rit; } virtual UInt getDim() { if (this->homogeneous) { - field_type_iterator tit = this->field.firstType( - this->spatial_dimension, this->ghost_type, this->element_kind); + auto tit = this->field.elementTypes( + this->spatial_dimension, this->ghost_type, this->element_kind).begin(); return this->getNbDataPerElem(*tit); } throw; return 0; } void setNbDataPerElem(const ElementTypeMap & nb_data) override { nb_data_per_elem = nb_data; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the ElementTypeMapArray embedded in the field const field_type & field; /// total number of elements UInt nb_total_element; /// the spatial dimension of the problem UInt spatial_dimension; /// whether this is a ghost field or not (for type selection) GhostType ghost_type; /// The element kind to operate on ElementKind element_kind; /// The number of data per element type ElementTypeMap nb_data_per_elem; }; /* -------------------------------------------------------------------------- */ #include "dumper_generic_elemental_field_tmpl.hh" /* -------------------------------------------------------------------------- */ __END_AKANTU_DUMPER__ } // akantu #endif /* __AKANTU_DUMPER_GENERIC_ELEMENTAL_FIELD_HH__ */ diff --git a/src/io/dumper/dumper_generic_elemental_field_tmpl.hh b/src/io/dumper/dumper_generic_elemental_field_tmpl.hh index f22b9be0f..6e00e2b66 100644 --- a/src/io/dumper/dumper_generic_elemental_field_tmpl.hh +++ b/src/io/dumper/dumper_generic_elemental_field_tmpl.hh @@ -1,69 +1,64 @@ /** * @file dumper_generic_elemental_field_tmpl.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Sep 02 2014 * @date last modification: Wed Nov 08 2017 * * @brief Implementation of the template functions of the ElementalField * * @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 . * */ /* -------------------------------------------------------------------------- */ template class iterator> void GenericElementalField::checkHomogeneity() { - using field_type_iterator = typename field_type::type_iterator; - using array_type = typename field_type::array_type; - - field_type_iterator tit; - field_type_iterator end; - - tit = field.firstType(spatial_dimension, ghost_type, element_kind); - end = field.lastType(spatial_dimension, ghost_type, element_kind); + auto types = field.elementTypes(spatial_dimension, ghost_type, element_kind); + auto tit = types.begin(); + auto end = types.end(); this->nb_total_element = 0; UInt nb_comp = 0; bool homogen = true; if (tit != end) { nb_comp = this->field(*tit, ghost_type).getNbComponent(); for (; tit != end; ++tit) { - const array_type & vect = this->field(*tit, ghost_type); - UInt nb_element = vect.size(); - UInt nb_comp_cur = vect.getNbComponent(); + const auto & vect = this->field(*tit, ghost_type); + auto nb_element = vect.size(); + auto nb_comp_cur = vect.getNbComponent(); if (homogen && nb_comp != nb_comp_cur) homogen = false; this->nb_total_element += nb_element; // this->nb_data_per_elem(*tit,this->ghost_type) = nb_comp_cur; } if (!homogen) nb_comp = 0; } this->homogeneous = homogen; } /* -------------------------------------------------------------------------- */ diff --git a/src/io/dumper/dumper_homogenizing_field.hh b/src/io/dumper/dumper_homogenizing_field.hh index c95bf75bb..daf4c9eb9 100644 --- a/src/io/dumper/dumper_homogenizing_field.hh +++ b/src/io/dumper/dumper_homogenizing_field.hh @@ -1,216 +1,216 @@ /** * @file dumper_homogenizing_field.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Sep 02 2014 * @date last modification: Wed Nov 08 2017 * * @brief description of field homogenizing field * * @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_DUMPER_HOMOGENIZING_FIELD_HH__ #define __AKANTU_DUMPER_HOMOGENIZING_FIELD_HH__ /* -------------------------------------------------------------------------- */ #include "dumper_compute.hh" /* -------------------------------------------------------------------------- */ namespace akantu { __BEGIN_AKANTU_DUMPER__ /* -------------------------------------------------------------------------- */ template inline type typeConverter(const type & input, __attribute__((unused)) Vector & res, __attribute__((unused)) UInt nb_data) { throw; return input; } /* -------------------------------------------------------------------------- */ template inline Matrix typeConverter(const Matrix & input, Vector & res, UInt nb_data) { Matrix tmp(res.storage(), input.rows(), nb_data / input.rows()); Matrix tmp2(tmp, true); return tmp2; } /* -------------------------------------------------------------------------- */ template -inline Vector typeConverter(__attribute__((unused)) - const Vector & input, - __attribute__((unused)) Vector & res, - __attribute__((unused)) UInt nb_data) { +inline Vector typeConverter(const Vector & , + Vector & res, + UInt) { return res; } /* -------------------------------------------------------------------------- */ template class AvgHomogenizingFunctor : public ComputeFunctor { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ using value_type = typename type::value_type; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: AvgHomogenizingFunctor(ElementTypeMap & nb_datas) { - ElementTypeMap::type_iterator tit = nb_datas.firstType(); - ElementTypeMap::type_iterator end = nb_datas.lastType(); + auto types = nb_datas.elementTypes(); + auto tit = types.begin(); + auto end = types.end(); nb_data = nb_datas(*tit); for (; tit != end; ++tit) if (nb_data != nb_datas(*tit)) throw; } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ type func(const type & d, Element /*global_index*/) override { Vector res(this->nb_data); if (d.size() % this->nb_data) throw; UInt nb_to_average = d.size() / this->nb_data; value_type * ptr = d.storage(); for (UInt i = 0; i < nb_to_average; ++i) { Vector tmp(ptr, this->nb_data); res += tmp; ptr += this->nb_data; } res /= nb_to_average; return typeConverter(d, res, this->nb_data); }; UInt getDim() override { return nb_data; }; UInt getNbComponent(UInt /*old_nb_comp*/) override { throw; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ /// The size of data: i.e. the size of the vector to be returned UInt nb_data; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class HomogenizerProxy { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: HomogenizerProxy() = default; public: inline static ComputeFunctorInterface * createHomogenizer(Field & field); template inline ComputeFunctorInterface * connectToField(T * field) { ElementTypeMap nb_components = field->getNbComponents(); using ret_type = typename T::types::return_type; return this->instantiateHomogenizer(nb_components); } template inline ComputeFunctorInterface * instantiateHomogenizer(ElementTypeMap & nb_components); }; /* -------------------------------------------------------------------------- */ template inline ComputeFunctorInterface * HomogenizerProxy::instantiateHomogenizer(ElementTypeMap & nb_components) { using Homogenizer = dumper::AvgHomogenizingFunctor; auto * foo = new Homogenizer(nb_components); return foo; } template <> inline ComputeFunctorInterface * HomogenizerProxy::instantiateHomogenizer>( __attribute__((unused)) ElementTypeMap & nb_components) { throw; return nullptr; } /* -------------------------------------------------------------------------- */ /// for connection to a FieldCompute template inline ComputeFunctorInterface * FieldCompute::connect(HomogenizerProxy & proxy) { return proxy.connectToField(this); } /* -------------------------------------------------------------------------- */ inline ComputeFunctorInterface * HomogenizerProxy::createHomogenizer(Field & field) { HomogenizerProxy homogenizer_proxy; return field.connect(homogenizer_proxy); } /* -------------------------------------------------------------------------- */ // inline ComputeFunctorInterface & createHomogenizer(Field & field){ // HomogenizerProxy::createHomogenizer(field); // throw; // ComputeFunctorInterface * ptr = NULL; // return *ptr; // } // /* -------------------------------------------------------------------------- // */ __END_AKANTU_DUMPER__ } // namespace akantu #endif /* __AKANTU_DUMPER_HOMOGENIZING_FIELD_HH__ */ diff --git a/src/io/mesh_io.cc b/src/io/mesh_io.cc index 23e96bb8a..e1f6a8214 100644 --- a/src/io/mesh_io.cc +++ b/src/io/mesh_io.cc @@ -1,147 +1,145 @@ /** * @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 if (ext == "inp") { t = _miot_abaqus; } 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(); case _miot_abaqus: 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()) { - for (Mesh::type_iterator type_it = mesh.firstType(); - type_it != mesh.lastType(); ++type_it) { - + for (auto type : mesh.elementTypes()) { auto & name_vec = - mesh.getDataPointer("physical_names", *type_it); + mesh.getDataPointer("physical_names", type); - const auto & tags_vec = mesh.getData(tag_name, *type_it); + 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); if (map_it == phys_name_map.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()) { stream << space << "Physical map:" << std::endl; for (auto & pair : phys_name_map) { stream << space << pair.first << ": " << pair.second << std::endl; } } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/io/mesh_io/mesh_io_diana.cc b/src/io/mesh_io/mesh_io_diana.cc index 625d0dfcd..2191b0a59 100644 --- a/src/io/mesh_io/mesh_io_diana.cc +++ b/src/io/mesh_io/mesh_io_diana.cc @@ -1,601 +1,597 @@ /** * @file mesh_io_diana.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * @author Alodie Schneuwly * * @date creation: Sat Mar 26 2011 * @date last modification: Tue Feb 20 2018 * * @brief handles diana meshes * * @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 #include /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "mesh_io_diana.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include namespace akantu { /* -------------------------------------------------------------------------- */ /* Methods Implentations */ /* -------------------------------------------------------------------------- */ MeshIODiana::MeshIODiana() { canReadSurface = true; canReadExtendedData = true; _diana_to_akantu_element_types["T9TM"] = _triangle_3; _diana_to_akantu_element_types["CT6CM"] = _triangle_6; _diana_to_akantu_element_types["Q12TM"] = _quadrangle_4; _diana_to_akantu_element_types["CQ8CM"] = _quadrangle_8; _diana_to_akantu_element_types["TP18L"] = _pentahedron_6; _diana_to_akantu_element_types["CTP45"] = _pentahedron_15; _diana_to_akantu_element_types["TE12L"] = _tetrahedron_4; _diana_to_akantu_element_types["HX24L"] = _hexahedron_8; _diana_to_akantu_element_types["CHX60"] = _hexahedron_20; _diana_to_akantu_mat_prop["YOUNG"] = "E"; _diana_to_akantu_mat_prop["DENSIT"] = "rho"; _diana_to_akantu_mat_prop["POISON"] = "nu"; std::map::iterator it; for (it = _diana_to_akantu_element_types.begin(); it != _diana_to_akantu_element_types.end(); ++it) { UInt nb_nodes = Mesh::getNbNodesPerElement(it->second); auto * tmp = new UInt[nb_nodes]; for (UInt i = 0; i < nb_nodes; ++i) { tmp[i] = i; } switch (it->second) { case _tetrahedron_10: tmp[8] = 9; tmp[9] = 8; break; case _pentahedron_15: tmp[0] = 2; tmp[1] = 8; tmp[2] = 0; tmp[3] = 6; tmp[4] = 1; tmp[5] = 7; tmp[6] = 11; tmp[7] = 9; tmp[8] = 10; tmp[9] = 5; tmp[10] = 14; tmp[11] = 3; tmp[12] = 12; tmp[13] = 4; tmp[14] = 13; break; case _hexahedron_20: tmp[0] = 5; tmp[1] = 16; tmp[2] = 4; tmp[3] = 19; tmp[4] = 7; tmp[5] = 18; tmp[6] = 6; tmp[7] = 17; tmp[8] = 13; tmp[9] = 12; tmp[10] = 15; tmp[11] = 14; tmp[12] = 1; tmp[13] = 8; tmp[14] = 0; tmp[15] = 11; tmp[16] = 3; tmp[17] = 10; tmp[18] = 2; tmp[19] = 9; break; default: // nothing to change break; } _read_order[it->second] = tmp; } } /* -------------------------------------------------------------------------- */ MeshIODiana::~MeshIODiana() = default; /* -------------------------------------------------------------------------- */ inline void my_getline(std::ifstream & infile, std::string & line) { std::getline(infile, line); // read the line size_t pos = line.find('\r'); /// remove the extra \r if needed line = line.substr(0, pos); } /* -------------------------------------------------------------------------- */ void MeshIODiana::read(const std::string & filename, Mesh & mesh) { AKANTU_DEBUG_IN(); MeshAccessor mesh_accessor(mesh); std::ifstream infile; infile.open(filename.c_str()); std::string line; UInt first_node_number = std::numeric_limits::max(); diana_element_number_to_elements.clear(); if (!infile.good()) { AKANTU_ERROR("Cannot open file " << filename); } while (infile.good()) { my_getline(infile, line); /// read all nodes if (line == "'COORDINATES'") { line = readCoordinates(infile, mesh, first_node_number); } /// read all elements if (line == "'ELEMENTS'") { line = readElements(infile, mesh, first_node_number); } /// read the material properties and write a .dat file if (line == "'MATERIALS'") { line = readMaterial(infile, filename); } /// read the material properties and write a .dat file if (line == "'GROUPS'") { line = readGroups(infile, mesh, first_node_number); } } infile.close(); mesh_accessor.setNbGlobalNodes(mesh.getNbNodes()); MeshUtils::fillElementToSubElementsData(mesh); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshIODiana::write(__attribute__((unused)) const std::string & filename, __attribute__((unused)) const Mesh & mesh) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readCoordinates(std::ifstream & infile, Mesh & mesh, UInt & first_node_number) { AKANTU_DEBUG_IN(); MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); std::string line; UInt index; Vector coord(3); do { my_getline(infile, line); if ("'ELEMENTS'" == line) break; std::stringstream sstr_node(line); sstr_node >> index >> coord(0) >> coord(1) >> coord(2); first_node_number = first_node_number < index ? first_node_number : index; nodes.push_back(coord); } while (true); AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ UInt MeshIODiana::readInterval(std::stringstream & line, std::set & interval) { UInt first; line >> first; if (line.fail()) { return 0; } interval.insert(first); UInt second; int dash; dash = line.get(); if (dash == '-') { line >> second; interval.insert(second); return 2; } if (line.fail()) line.clear(std::ios::eofbit); // in case of get at end of the line else line.unget(); return 1; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readGroups(std::ifstream & infile, Mesh & mesh, UInt first_node_number) { AKANTU_DEBUG_IN(); std::string line; my_getline(infile, line); bool reading_nodes_group = false; while (line != "'SUPPORTS'") { if (line == "NODES") { reading_nodes_group = true; my_getline(infile, line); } if (line == "ELEMEN") { reading_nodes_group = false; my_getline(infile, line); } auto * str = new std::stringstream(line); UInt id; std::string name; char c; *str >> id >> name >> c; auto * list_ids = new Array(0, 1, name); UInt s = 1; bool end = false; while (!end) { while (!str->eof() && s != 0) { std::set interval; s = readInterval(*str, interval); auto it = interval.begin(); if (s == 1) list_ids->push_back(*it); if (s == 2) { UInt first = *it; ++it; UInt second = *it; for (UInt i = first; i <= second; ++i) { list_ids->push_back(i); } } } if (str->fail()) end = true; else { my_getline(infile, line); delete str; str = new std::stringstream(line); } } delete str; if (reading_nodes_group) { NodeGroup & ng = mesh.createNodeGroup(name); for (UInt i = 0; i < list_ids->size(); ++i) { UInt node = (*list_ids)(i)-first_node_number; ng.add(node, false); } delete list_ids; } else { ElementGroup & eg = mesh.createElementGroup(name); for (UInt i = 0; i < list_ids->size(); ++i) { Element & elem = diana_element_number_to_elements[(*list_ids)(i)]; if (elem.type != _not_defined) eg.add(elem, false, false); } eg.optimize(); delete list_ids; } my_getline(infile, line); } AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readElements(std::ifstream & infile, Mesh & mesh, UInt first_node_number) { AKANTU_DEBUG_IN(); std::string line; my_getline(infile, line); if ("CONNECTIVITY" == line) { line = readConnectivity(infile, mesh, first_node_number); } /// read the line corresponding to the materials if ("MATERIALS" == line) { line = readMaterialElement(infile, mesh); } AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readConnectivity(std::ifstream & infile, Mesh & mesh, UInt first_node_number) { AKANTU_DEBUG_IN(); MeshAccessor mesh_accessor(mesh); Int index; std::string lline; std::string diana_type; ElementType akantu_type, akantu_type_old = _not_defined; Array * connectivity = nullptr; UInt node_per_element = 0; Element elem; UInt * read_order = nullptr; while (true) { my_getline(infile, lline); // std::cerr << lline << std::endl; std::stringstream sstr_elem(lline); if (lline == "MATERIALS") break; /// traiter les coordonnees sstr_elem >> index; sstr_elem >> diana_type; akantu_type = _diana_to_akantu_element_types[diana_type]; if (akantu_type == _not_defined) continue; if (akantu_type != akantu_type_old) { connectivity = &(mesh_accessor.getConnectivity(akantu_type)); node_per_element = connectivity->getNbComponent(); akantu_type_old = akantu_type; read_order = _read_order[akantu_type]; } Vector local_connect(node_per_element); // used if element is written on two lines UInt j_last = 0; for (UInt j = 0; j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; // check s'il y a pas plus rien après un certain point if (sstr_elem.fail()) { sstr_elem.clear(); sstr_elem.ignore(); break; } node_index -= first_node_number; local_connect(read_order[j]) = node_index; j_last = j; } // check if element is written in two lines if (j_last != (node_per_element - 1)) { // if this is the case, read on more line my_getline(infile, lline); std::stringstream sstr_elem(lline); for (UInt j = (j_last + 1); j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; node_index -= first_node_number; local_connect(read_order[j]) = node_index; } } connectivity->push_back(local_connect); elem.type = akantu_type; elem.element = connectivity->size() - 1; diana_element_number_to_elements[index] = elem; akantu_number_to_diana_number[elem] = index; } AKANTU_DEBUG_OUT(); return lline; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readMaterialElement(std::ifstream & infile, Mesh & mesh) { AKANTU_DEBUG_IN(); std::string line; - std::stringstream sstr_tag_name; - sstr_tag_name << "tag_" << 0; - - Mesh::type_iterator it = mesh.firstType(); - Mesh::type_iterator end = mesh.lastType(); - for (; it != end; ++it) { - UInt nb_element = mesh.getNbElement(*it); - mesh.getDataPointer("material", *it, _not_ghost, 1) + + for (auto type : mesh.elementTypes()) { + UInt nb_element = mesh.getNbElement(type); + mesh.getDataPointer("material", type, _not_ghost, 1) .resize(nb_element); } my_getline(infile, line); while (line != "'MATERIALS'") { line = line.substr(line.find('/') + 1, std::string::npos); // erase the first slash / of the line char tutu[250] = {'\0'}; strncpy(tutu, line.c_str(), 249); AKANTU_DEBUG_WARNING("reading line " << line); Array temp_id(0, 2); UInt mat; while (true) { std::stringstream sstr_intervals_elements(line); Vector id(2); char temp; while (sstr_intervals_elements.good()) { sstr_intervals_elements >> id(0) >> temp >> id(1); // >> "/" >> mat; if (!sstr_intervals_elements.fail()) temp_id.push_back(id); } if (sstr_intervals_elements.fail()) { sstr_intervals_elements.clear(); sstr_intervals_elements.ignore(); sstr_intervals_elements >> mat; break; } my_getline(infile, line); } // loop over elements // UInt * temp_id_val = temp_id.storage(); for (UInt i = 0; i < temp_id.size(); ++i) for (UInt j = temp_id(i, 0); j <= temp_id(i, 1); ++j) { Element & element = diana_element_number_to_elements[j]; if (element.type == _not_defined) continue; UInt elem = element.element; ElementType type = element.type; Array & data = mesh.getDataPointer("material", type, _not_ghost); data(elem) = mat; } my_getline(infile, line); } AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readMaterial(std::ifstream & infile, const std::string & filename) { AKANTU_DEBUG_IN(); std::stringstream mat_file_name; mat_file_name << "material_" << filename; std::ofstream material_file; material_file.open(mat_file_name.str().c_str()); // mat_file_name.str()); UInt mat_index; std::string line; bool first_mat = true; bool end = false; UInt mat_id = 0; using MatProp = std::map; MatProp mat_prop; do { my_getline(infile, line); std::stringstream sstr_material(line); if (("'GROUPS'" == line) || ("'END'" == line)) { if (!mat_prop.empty()) { material_file << "material elastic [" << std::endl; material_file << "\tname = material" << ++mat_id << std::endl; for (auto it = mat_prop.begin(); it != mat_prop.end(); ++it) material_file << "\t" << it->first << " = " << it->second << std::endl; material_file << "]" << std::endl; mat_prop.clear(); } end = true; } else { /// traiter les caractéristiques des matériaux sstr_material >> mat_index; if (!sstr_material.fail()) { if (!first_mat) { if (!mat_prop.empty()) { material_file << "material elastic [" << std::endl; material_file << "\tname = material" << ++mat_id << std::endl; for (auto it = mat_prop.begin(); it != mat_prop.end(); ++it) material_file << "\t" << it->first << " = " << it->second << std::endl; material_file << "]" << std::endl; mat_prop.clear(); } } first_mat = false; } else { sstr_material.clear(); } std::string prop_name; sstr_material >> prop_name; std::map::iterator it; it = _diana_to_akantu_mat_prop.find(prop_name); if (it != _diana_to_akantu_mat_prop.end()) { Real value; sstr_material >> value; mat_prop[it->second] = value; } else { AKANTU_DEBUG_INFO("In material reader, property " << prop_name << "not recognized"); } } } while (!end); AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/io/parser/parser.hh b/src/io/parser/parser.hh index ea07a640b..b724fc7ea 100644 --- a/src/io/parser/parser.hh +++ b/src/io/parser/parser.hh @@ -1,509 +1,509 @@ /** * @file parser.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Fri Dec 08 2017 * * @brief File parser interface * * @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_random_generator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PARSER_HH__ #define __AKANTU_PARSER_HH__ namespace akantu { #ifndef SWIG // clang-format off #define AKANTU_SECTION_TYPES \ (cohesive_inserter) \ (contact) \ (embedded_interface) \ (friction) \ (global) \ (heat) \ (integration_scheme) \ (material) \ (mesh) \ (model) \ (model_solver) \ (neighborhood) \ (neighborhoods) \ (non_linear_solver) \ (non_local) \ (rules) \ (solver) \ (time_step_solver) \ (user) \ (weight_function) \ (not_defined) // clang-format on /// Defines the possible section types -AKANTU_ENUM_DECLARE(ParserType, AKANTU_SECTION_TYPES) -AKANTU_ENUM_OUTPUT_STREAM(ParserType, AKANTU_SECTION_TYPES) -AKANTU_ENUM_INPUT_STREAM(ParserType, AKANTU_SECTION_TYPES) +AKANTU_CLASS_ENUM_DECLARE(ParserType, AKANTU_SECTION_TYPES) +AKANTU_CLASS_ENUM_OUTPUT_STREAM(ParserType, AKANTU_SECTION_TYPES) +AKANTU_CLASS_ENUM_INPUT_STREAM(ParserType, AKANTU_SECTION_TYPES) /// Defines the possible search contexts/scopes (for parameter search) enum ParserParameterSearchCxt { _ppsc_current_scope = 0x1, _ppsc_parent_scope = 0x2, _ppsc_current_and_parent_scope = 0x3 }; #endif /* ------------------------------------------------------------------------ */ /* Parameters Class */ /* ------------------------------------------------------------------------ */ class ParserSection; /// @brief The ParserParameter objects represent the end of tree branches as /// they /// are the different informations contained in the input file. class ParserParameter { public: ParserParameter() : name(std::string()), value(std::string()), dbg_filename(std::string()) { } ParserParameter(const std::string & name, const std::string & value, const ParserSection & parent_section) : parent_section(&parent_section), name(name), value(value), dbg_filename(std::string()) {} ParserParameter(const ParserParameter & param) = default; virtual ~ParserParameter() = default; /// Get parameter name const std::string & getName() const { return name; } /// Get parameter value const std::string & getValue() const { return value; } /// Set info for debug output void setDebugInfo(const std::string & filename, UInt line, UInt column) { dbg_filename = filename; dbg_line = line; dbg_column = column; } template inline operator T() const; // template inline operator Vector() const; // template inline operator Matrix() const; /// Print parameter info in stream void printself(std::ostream & stream, __attribute__((unused)) unsigned int indent = 0) const { stream << name << ": " << value << " (" << dbg_filename << ":" << dbg_line << ":" << dbg_column << ")"; } private: void setParent(const ParserSection & sect) { parent_section = § } friend class ParserSection; private: /// Pointer to the parent section const ParserSection * parent_section{nullptr}; /// Name of the parameter std::string name; /// Value of the parameter std::string value; /// File for debug output std::string dbg_filename; /// Position of parameter in parsed file UInt dbg_line, dbg_column; }; /* ------------------------------------------------------------------------ */ /* Sections Class */ /* ------------------------------------------------------------------------ */ /// ParserSection represents a branch of the parsing tree. class ParserSection { public: using SubSections = std::multimap; using Parameters = std::map; private: using const_section_iterator_ = SubSections::const_iterator; public: /* ------------------------------------------------------------------------ */ /* SubSection iterator */ /* ------------------------------------------------------------------------ */ /// Iterator on sections class const_section_iterator { public: using iterator_category = std::forward_iterator_tag; using value_type = ParserSection; using pointer = ParserSection *; using reference = ParserSection &; const_section_iterator() = default; const_section_iterator(const const_section_iterator_ & it) : it(it) {} const_section_iterator(const const_section_iterator & other) = default; const_section_iterator & operator=(const const_section_iterator & other) = default; const ParserSection & operator*() const { return it->second; } const ParserSection * operator->() const { return &(it->second); } bool operator==(const const_section_iterator & other) const { return it == other.it; } bool operator!=(const const_section_iterator & other) const { return it != other.it; } const_section_iterator & operator++() { ++it; return *this; } const_section_iterator operator++(int) { const_section_iterator tmp = *this; operator++(); return tmp; } private: const_section_iterator_ it; }; /* ------------------------------------------------------------------------ */ /* Parameters iterator */ /* ------------------------------------------------------------------------ */ /// Iterator on parameters class const_parameter_iterator { public: const_parameter_iterator(const const_parameter_iterator & other) = default; const_parameter_iterator(const Parameters::const_iterator & it) : it(it) {} const_parameter_iterator & operator=(const const_parameter_iterator & other) { if (this != &other) { it = other.it; } return *this; } const ParserParameter & operator*() const { return it->second; } const ParserParameter * operator->() { return &(it->second); }; bool operator==(const const_parameter_iterator & other) const { return it == other.it; } bool operator!=(const const_parameter_iterator & other) const { return it != other.it; } const_parameter_iterator & operator++() { ++it; return *this; } const_parameter_iterator operator++(int) { const_parameter_iterator tmp = *this; operator++(); return tmp; } private: Parameters::const_iterator it; }; /* ---------------------------------------------------------------------- */ ParserSection() : name(std::string()) {} ParserSection(const std::string & name, ParserType type) : parent_section(nullptr), name(name), type(type) {} ParserSection(const std::string & name, ParserType type, const std::string & option, const ParserSection & parent_section) : parent_section(&parent_section), name(name), type(type), option(option) {} ParserSection(const ParserSection & section) : parent_section(section.parent_section), name(section.name), type(section.type), option(section.option), parameters(section.parameters), sub_sections_by_type(section.sub_sections_by_type) { setChldrenPointers(); } ParserSection & operator=(const ParserSection & other) { if (&other != this) { parent_section = other.parent_section; name = other.name; type = other.type; option = other.option; parameters = other.parameters; sub_sections_by_type = other.sub_sections_by_type; setChldrenPointers(); } return *this; } virtual ~ParserSection(); virtual void printself(std::ostream & stream, unsigned int indent = 0) const; /* ---------------------------------------------------------------------- */ /* Creation functions */ /* ---------------------------------------------------------------------- */ public: ParserParameter & addParameter(const ParserParameter & param); ParserSection & addSubSection(const ParserSection & section); protected: /// Clean ParserSection content void clean() { parameters.clear(); sub_sections_by_type.clear(); } private: void setChldrenPointers() { for (auto && param_pair : this->parameters) param_pair.second.setParent(*this); for (auto && sub_sect_pair : this->sub_sections_by_type) sub_sect_pair.second.setParent(*this); } /* ---------------------------------------------------------------------- */ /* Accessors */ /* ---------------------------------------------------------------------- */ public: #ifndef SWIG class SubSectionsRange : public std::pair { public: SubSectionsRange(const const_section_iterator & first, const const_section_iterator & second) : std::pair(first, second) {} auto begin() { return this->first; } auto end() { return this->second; } }; /// Get begin and end iterators on subsections of certain type auto getSubSections(ParserType type = ParserType::_not_defined) const { if (type != ParserType::_not_defined) { auto range = sub_sections_by_type.equal_range(type); return SubSectionsRange(range.first, range.second); } else { return SubSectionsRange(sub_sections_by_type.begin(), sub_sections_by_type.end()); } } /// Get number of subsections of certain type UInt getNbSubSections(ParserType type = ParserType::_not_defined) const { if (type != ParserType::_not_defined) { return this->sub_sections_by_type.count(type); } else { return this->sub_sections_by_type.size(); } } /// Get begin and end iterators on parameters auto getParameters() const { return std::pair( parameters.begin(), parameters.end()); } #endif /* ---------------------------------------------------------------------- */ /// Get parameter within specified context const ParserParameter & getParameter( const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if (search_ctx & _ppsc_current_scope) it = parameters.find(name); if (it == parameters.end()) { if ((search_ctx & _ppsc_parent_scope) && parent_section) return parent_section->getParameter(name, search_ctx); else { AKANTU_SILENT_EXCEPTION( "The parameter " << name << " has not been found in the specified context"); } } return it->second; } /* ------------------------------------------------------------------------ */ /// Get parameter within specified context, with a default value in case the /// parameter does not exists template const T getParameter( const std::string & name, const T & default_value, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { try { T tmp = this->getParameter(name, search_ctx); return tmp; } catch (debug::Exception &) { return default_value; } } /* ------------------------------------------------------------------------ */ /// Check if parameter exists within specified context bool hasParameter( const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if (search_ctx & _ppsc_current_scope) it = parameters.find(name); if (it == parameters.end()) { if ((search_ctx & _ppsc_parent_scope) && parent_section) return parent_section->hasParameter(name, search_ctx); else { return false; } } return true; } /* -------------------------------------------------------------------------- */ /// Get value of given parameter in context template T getParameterValue( const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { const ParserParameter & tmp_param = getParameter(name, search_ctx); T t = tmp_param; return t; } /* -------------------------------------------------------------------------- */ /// Get section name const std::string getName() const { return name; } /// Get section type ParserType getType() const { return type; } /// Get section option const std::string getOption(const std::string & def = "") const { return option != "" ? option : def; } protected: void setParent(const ParserSection & sect) { parent_section = § } /* ---------------------------------------------------------------------- */ /* Members */ /* ---------------------------------------------------------------------- */ private: /// Pointer to the parent section const ParserSection * parent_section{nullptr}; /// Name of section std::string name; /// Type of section, see AKANTU_SECTION_TYPES ParserType type{ParserType::_not_defined}; /// Section option std::string option; /// Map of parameters in section Parameters parameters; /// Multi-map of subsections SubSections sub_sections_by_type; }; /* ------------------------------------------------------------------------ */ /* Parser Class */ /* ------------------------------------------------------------------------ */ /// Root of parsing tree, represents the global ParserSection class Parser : public ParserSection { public: Parser() : ParserSection("global", ParserType::_global) {} void parse(const std::string & filename); std::string getLastParsedFile() const; static bool isPermissive() { return permissive_parser; } public: /// Parse real scalar static Real parseReal(const std::string & value, const ParserSection & section); /// Parse real vector static Vector parseVector(const std::string & value, const ParserSection & section); /// Parse real matrix static Matrix parseMatrix(const std::string & value, const ParserSection & section); #ifndef SWIG /// Parse real random parameter static RandomParameter parseRandomParameter(const std::string & value, const ParserSection & section); #endif protected: /// General parse function template static T parseType(const std::string & value, Grammar & grammar); protected: // friend class Parsable; static bool permissive_parser; std::string last_parsed_file; }; inline std::ostream & operator<<(std::ostream & stream, const ParserParameter & _this) { _this.printself(stream); return stream; } inline std::ostream & operator<<(std::ostream & stream, const ParserSection & section) { section.printself(stream); return stream; } } // akantu namespace std { template <> struct iterator_traits<::akantu::Parser::const_section_iterator> { using iterator_category = input_iterator_tag; using value_type = ::akantu::ParserParameter; using difference_type = ptrdiff_t; using pointer = const ::akantu::ParserParameter *; using reference = const ::akantu::ParserParameter &; }; } #include "parser_tmpl.hh" #endif /* __AKANTU_PARSER_HH__ */ diff --git a/src/mesh/element_group.cc b/src/mesh/element_group.cc index ded48bb53..c29d194a9 100644 --- a/src/mesh/element_group.cc +++ b/src/mesh/element_group.cc @@ -1,203 +1,188 @@ /** * @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); #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 (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - - GhostType ghost_type = *gt; - - type_iterator it = - other_group.firstType(_all_dimensions, ghost_type, _ek_not_defined); - type_iterator last = - other_group.lastType(_all_dimensions, ghost_type, _ek_not_defined); - - for (; it != last; ++it) { - ElementType type = *it; + 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 (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - GhostType ghost_type = *gt; - ElementList::type_iterator it = - elements.firstType(_all_dimensions, ghost_type); - ElementList::type_iterator last = - elements.lastType(_all_dimensions, ghost_type); - for (; it != last; ++it) { - Array & els = elements(*it, ghost_type); + 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(); } /* -------------------------------------------------------------------------- */ -} // akantu +} // namespace akantu diff --git a/src/mesh/element_group.hh b/src/mesh/element_group.hh index 64615adaf..653326c28 100644 --- a/src/mesh/element_group.hh +++ b/src/mesh/element_group.hh @@ -1,190 +1,193 @@ /** * @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 */ /* ------------------------------------------------------------------------ */ 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; #ifndef SWIG 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(); 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(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/element_group_inline_impl.cc b/src/mesh/element_group_inline_impl.cc index 7ae8b3750..adf88a2c6 100644 --- a/src/mesh/element_group_inline_impl.cc +++ b/src/mesh/element_group_inline_impl.cc @@ -1,146 +1,146 @@ /** * @file element_group_inline_impl.cc * * @author Dana Christen * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Sun Aug 13 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 "element_group.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_GROUP_INLINE_IMPL_CC__ #define __AKANTU_ELEMENT_GROUP_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline void ElementGroup::add(const Element & el, bool add_nodes, bool check_for_duplicate) { this->add(el.type, el.element, el.ghost_type, add_nodes, check_for_duplicate); } /* -------------------------------------------------------------------------- */ inline void ElementGroup::add(const ElementType & type, UInt element, const GhostType & ghost_type, bool add_nodes, bool check_for_duplicate) { addElement(type, element, ghost_type); if (add_nodes) { Array::const_vector_iterator it = mesh.getConnectivity(type, ghost_type) .begin(mesh.getNbNodesPerElement(type)) + element; const Vector & conn = *it; for (UInt i = 0; i < conn.size(); ++i) addNode(conn[i], check_for_duplicate); } } /* -------------------------------------------------------------------------- */ inline void ElementGroup::addNode(UInt node_id, bool check_for_duplicate) { node_group.add(node_id, check_for_duplicate); } /* -------------------------------------------------------------------------- */ inline void ElementGroup::removeNode(UInt node_id) { node_group.remove(node_id); } /* -------------------------------------------------------------------------- */ inline void ElementGroup::addElement(const ElementType & elem_type, UInt elem_id, const GhostType & ghost_type) { if (!(elements.exists(elem_type, ghost_type))) { elements.alloc(0, 1, elem_type, ghost_type); } elements(elem_type, ghost_type).push_back(elem_id); this->dimension = UInt( std::max(Int(this->dimension), Int(mesh.getSpatialDimension(elem_type)))); } /* -------------------------------------------------------------------------- */ inline UInt ElementGroup::getNbNodes() const { return node_group.size(); } /* -------------------------------------------------------------------------- */ inline ElementGroup::type_iterator ElementGroup::firstType(UInt dim, const GhostType & ghost_type, const ElementKind & kind) const { - return elements.firstType(dim, ghost_type, kind); + return elements.elementTypes(dim, ghost_type, kind).begin(); } /* -------------------------------------------------------------------------- */ inline ElementGroup::type_iterator ElementGroup::lastType(UInt dim, const GhostType & ghost_type, const ElementKind & kind) const { - return elements.lastType(dim, ghost_type, kind); + return elements.elementTypes(dim, ghost_type, kind).end(); } /* -------------------------------------------------------------------------- */ inline ElementGroup::const_element_iterator ElementGroup::begin(const ElementType & type, const GhostType & ghost_type) const { if (elements.exists(type, ghost_type)) { return elements(type, ghost_type).begin(); } else { return empty_elements.begin(); } } /* -------------------------------------------------------------------------- */ inline ElementGroup::const_element_iterator ElementGroup::end(const ElementType & type, const GhostType & ghost_type) const { if (elements.exists(type, ghost_type)) { return elements(type, ghost_type).end(); } else { return empty_elements.end(); } } /* -------------------------------------------------------------------------- */ inline const Array & ElementGroup::getElements(const ElementType & type, const GhostType & ghost_type) const { if (elements.exists(type, ghost_type)) { return elements(type, ghost_type); } else { return empty_elements; } } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_ELEMENT_GROUP_INLINE_IMPL_CC__ */ diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh index 4853569a7..cbf2c37cf 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,449 +1,447 @@ /** * @file element_type_map.hh * * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Tue Feb 20 2018 * * @brief storage class by element type * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_memory.hh" #include "aka_named_argument.hh" #include "element.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_TYPE_MAP_HH__ #define __AKANTU_ELEMENT_TYPE_MAP_HH__ namespace akantu { class FEEngine; } // namespace akantu namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(all_ghost_types); DECLARE_NAMED_ARGUMENT(default_value); DECLARE_NAMED_ARGUMENT(element_kind); DECLARE_NAMED_ARGUMENT(ghost_type); DECLARE_NAMED_ARGUMENT(nb_component); DECLARE_NAMED_ARGUMENT(nb_component_functor); DECLARE_NAMED_ARGUMENT(with_nb_element); DECLARE_NAMED_ARGUMENT(with_nb_nodes_per_element); DECLARE_NAMED_ARGUMENT(spatial_dimension); DECLARE_NAMED_ARGUMENT(do_not_default); } // namespace template class ElementTypeMap; /* -------------------------------------------------------------------------- */ /* ElementTypeMapBase */ /* -------------------------------------------------------------------------- */ /// Common non templated base class for the ElementTypeMap class class ElementTypeMapBase { public: virtual ~ElementTypeMapBase() = default; }; /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template class ElementTypeMap : public ElementTypeMapBase { public: ElementTypeMap(); ~ElementTypeMap() override; inline static std::string printType(const SupportType & type, const GhostType & ghost_type); /*! Tests whether a type is present in the object * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return true if the type is present. */ inline bool exists(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline const Stored & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline Stored & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. THIS METHOD IS * NOT ARRAY SAFE, when using ElementTypeMapArray, use setArray instead * @param data to insert * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ template inline Stored & operator()(U && insertee, const SupportType & type, const GhostType & ghost_type = _not_ghost); public: /// print helper virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ /*! iterator allows to iterate over type-data pairs of the map. The interface * expects the SupportType to be ElementType. */ using DataMap = std::map; + + /// helper class to use in range for constructions class type_iterator : private std::iterator { public: using value_type = const SupportType; using pointer = const SupportType *; using reference = const SupportType &; protected: using DataMapIterator = typename ElementTypeMap::DataMap::const_iterator; public: type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek); type_iterator(const type_iterator & it); type_iterator() = default; inline reference operator*(); inline reference operator*() const; inline type_iterator & operator++(); type_iterator operator++(int); inline bool operator==(const type_iterator & other) const; inline bool operator!=(const type_iterator & other) const; type_iterator & operator=(const type_iterator & other); private: DataMapIterator list_begin; DataMapIterator list_end; UInt dim; ElementKind kind; }; /// helper class to use in range for constructions class ElementTypesIteratorHelper { public: using Container = ElementTypeMap; using iterator = typename Container::type_iterator; ElementTypesIteratorHelper(const Container & container, UInt dim, GhostType ghost_type, ElementKind kind) : container(std::cref(container)), dim(dim), ghost_type(ghost_type), kind(kind) {} template ElementTypesIteratorHelper(const Container & container, use_named_args_t, pack &&... _pack) : ElementTypesIteratorHelper( container, OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), OPTIONAL_NAMED_ARG(ghost_type, _not_ghost), OPTIONAL_NAMED_ARG(element_kind, _ek_regular)) {} ElementTypesIteratorHelper(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(ElementTypesIteratorHelper &&) = default; - iterator begin() { - return container.get().firstType(dim, ghost_type, kind); - } - iterator end() { return container.get().lastType(dim, ghost_type, kind); } + iterator begin(); + iterator end(); private: std::reference_wrapper container; UInt dim; GhostType ghost_type; ElementKind kind; }; private: ElementTypesIteratorHelper elementTypesImpl(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const; template ElementTypesIteratorHelper elementTypesImpl(const use_named_args_t & /*unused*/, pack &&... _pack) const; public: /*! * \param _spatial_dimension optional: filter for elements of given spatial * dimension * \param _ghost_type optional: filter for a certain ghost_type * \param _element_kind optional: filter for elements of given kind */ template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(use_named_args, std::forward(_pack)...); } template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(std::forward(_pack)...); } -public: /*! Get an iterator to the beginning of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the first stored data matching the filters * or an iterator to the end of the map if none match*/ - inline type_iterator firstType(UInt dim = _all_dimensions, - GhostType ghost_type = _not_ghost, - ElementKind kind = _ek_not_defined) const; + [[deprecated("Use elementTypes instead")]] inline type_iterator + firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, + ElementKind kind = _ek_not_defined) const; /*! Get an iterator to the end of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the last stored data matching the filters * or an iterator to the end of the map if none match */ - inline type_iterator lastType(UInt dim = _all_dimensions, - GhostType ghost_type = _not_ghost, - ElementKind kind = _ek_not_defined) const; + [[deprecated("Use elementTypes instead")]] inline type_iterator + lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, + ElementKind kind = _ek_not_defined) const; -protected: /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline DataMap & getData(GhostType ghost_type); /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline const DataMap & getData(GhostType ghost_type) const; /* ------------------------------------------------------------------------ */ protected: DataMap data; DataMap ghost_data; }; /* -------------------------------------------------------------------------- */ /* Some typedefs */ /* -------------------------------------------------------------------------- */ template class ElementTypeMapArray : public ElementTypeMap *, SupportType>, public Memory { public: using type = T; using array_type = Array; protected: using parent = ElementTypeMap *, SupportType>; using DataMap = typename parent::DataMap; public: using type_iterator = typename parent::type_iterator; /// standard assigment (copy) operator void operator=(const ElementTypeMapArray &) = delete; ElementTypeMapArray(const ElementTypeMapArray &) = delete; /// explicit copy void copy(const ElementTypeMapArray & other); /*! Constructor * @param id optional: identifier (string) * @param parent_id optional: parent identifier. for organizational purposes * only * @param memory_id optional: choose a specific memory, defaults to memory 0 */ ElementTypeMapArray(const ID & id = "by_element_type_array", const ID & parent_id = "no_parent", const MemoryID & memory_id = 0) : parent(), Memory(parent_id + ":" + id, memory_id), name(id){}; /*! allocate memory for a new array * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map * @param ghost_type whether to add the field to the data map or the * ghost_data map * @return a reference to the allocated array */ inline Array & alloc(UInt size, UInt nb_component, const SupportType & type, const GhostType & ghost_type, const T & default_value = T()); /*! allocate memory for a new array in both the data and the ghost_data map * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map*/ inline void alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value = T()); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a reference to the array */ inline const Array & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /// access the data of an element, this combine the map and array accessor inline const T & operator()(const Element & element, UInt component = 0) const; /// access the data of an element, this combine the map and array accessor inline T & operator()(const Element & element, UInt component = 0); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a const reference to the array */ inline Array & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @param vect the vector to include into the map * @return stored data corresponding to type. */ inline void setArray(const SupportType & type, const GhostType & ghost_type, const Array & vect); /*! frees all memory related to the data*/ inline void free(); /*! set all values in the ElementTypeMap to zero*/ inline void clear(); /*! set all values in the ElementTypeMap to value */ template inline void set(const ST & value); /*! deletes and reorders entries in the stored arrays * @param new_numbering a ElementTypeMapArray of new indices. UInt(-1) * indicates * deleted entries. */ inline void onElementsRemoved(const ElementTypeMapArray & new_numbering); /// text output helper void printself(std::ostream & stream, int indent = 0) const override; /*! set the id * @param id the new name */ inline void setID(const ID & id) { this->id = id; } ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const { ElementTypeMap nb_components; for (auto & type : this->elementTypes(dim, ghost_type, kind)) { UInt nb_comp = (*this)(type, ghost_type).getNbComponent(); nb_components(type, ghost_type) = nb_comp; } return nb_components; } /* ------------------------------------------------------------------------ */ /* more evolved allocators */ /* ------------------------------------------------------------------------ */ public: /// initialize the arrays in accordance to a functor template void initialize(const Func & f, const T & default_value, bool do_not_default); /// initialize with sizes and number of components in accordance of a mesh /// content template void initialize(const Mesh & mesh, pack &&... _pack); /// initialize with sizes and number of components in accordance of a fe /// engine content (aka integration points) template void initialize(const FEEngine & fe_engine, pack &&... _pack); /* ------------------------------------------------------------------------ */ /* Accesssors */ /* ------------------------------------------------------------------------ */ public: /// get the name of the internal field AKANTU_GET_MACRO(Name, name, ID); /// name of the elment type map: e.g. connectivity, grad_u ID name; }; /// to store data Array by element type using ElementTypeMapReal = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapInt = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapUInt = ElementTypeMapArray; /// Map of data of type UInt stored in a mesh using UIntDataMap = std::map *>; using ElementTypeMapUIntDataMap = ElementTypeMap; } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_HH__ */ diff --git a/src/mesh/element_type_map_filter.hh b/src/mesh/element_type_map_filter.hh index 834469672..10cf54011 100644 --- a/src/mesh/element_type_map_filter.hh +++ b/src/mesh/element_type_map_filter.hh @@ -1,307 +1,300 @@ /** * @file element_type_map_filter.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Sep 02 2014 * @date last modification: Sun Dec 03 2017 * * @brief Filtered version based on a an akantu::ElementGroup of a * akantu::ElementTypeMap * * @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_BY_ELEMENT_TYPE_FILTER_HH__ #define __AKANTU_BY_ELEMENT_TYPE_FILTER_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* ArrayFilter */ /* -------------------------------------------------------------------------- */ template class ArrayFilter { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ public: /// standard iterator template class iterator { inline bool operator!=(__attribute__((unused)) iterator & other) { throw; }; inline bool operator==(__attribute__((unused)) iterator & other) { throw; }; inline iterator & operator++() { throw; }; inline T operator*() { throw; return T(); }; }; /// const iterator template