diff --git a/examples/io/dumper/dumpable_interface.cc b/examples/io/dumper/dumpable_interface.cc index dc4edaa50..fe3505088 100644 --- a/examples/io/dumper/dumpable_interface.cc +++ b/examples/io/dumper/dumpable_interface.cc @@ -1,185 +1,186 @@ /** * @file dumpable_interface.cc * * @author Fabian Barras * @author Nicolas Richart * * @date creation: Mon Aug 17 2015 * @date last modification: Mon Aug 31 2015 * * @brief Example of dumper::Dumpable interface. * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "group_manager_inline_impl.cc" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include "locomotive_tools.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { /* In this example, we present dumper::Dumpable which is an interface for other classes who want to dump themselves. Several classes of Akantu inheritate from Dumpable (Model, Mesh, ...). In this example we reproduce the same tasks as example_dumper_low_level.cc using this time Dumpable interface inherted by Mesh, NodeGroup and ElementGroup. It is then advised to read first example_dumper_low_level.cc. */ initialize(argc, argv); // To start let us load the swiss train mesh and its mesh data information. UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); mesh.read("swiss_train.msh"); /* swiss_train.msh has the following physical groups that can be viewed with GMSH: "$MeshFormat 2.2 0 8 $EndMeshFormat $PhysicalNames 6 2 1 "red" 2 2 "white" 2 3 "lwheel_1" 2 4 "lwheel_2" 2 5 "rwheel_2" 2 6 "rwheel_1" $EndPhysicalNames ..." */ // Grouping nodes and elements belonging to train wheels (=four mesh data). ElementGroup & wheels_elements = mesh.createElementGroup("wheels", spatial_dimension); wheels_elements.append(mesh.getElementGroup("lwheel_1")); wheels_elements.append(mesh.getElementGroup("lwheel_2")); wheels_elements.append(mesh.getElementGroup("rwheel_1")); wheels_elements.append(mesh.getElementGroup("rwheel_2")); const Array & lnode_1 = (mesh.getElementGroup("lwheel_1")).getNodes(); const Array & lnode_2 = (mesh.getElementGroup("lwheel_2")).getNodes(); const Array & rnode_1 = (mesh.getElementGroup("rwheel_1")).getNodes(); const Array & rnode_2 = (mesh.getElementGroup("rwheel_2")).getNodes(); Array & node = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); // This time a 2D Array is created and a padding size of 3 is passed to // NodalField in order to warp train deformation on Paraview. Array displacement(nb_nodes, spatial_dimension); // Create an ElementTypeMapArray for the colour ElementTypeMapArray colour("colour"); colour.initialize(mesh, _with_nb_element = true); /* ------------------------------------------------------------------------ */ /* Creating dumpers */ /* ------------------------------------------------------------------------ */ // Create dumper for the complete mesh and register it as default dumper. DumperParaview dumper("train", "./paraview/dumpable", false); mesh.registerExternalDumper(dumper, "train", true); mesh.addDumpMesh(mesh); // The dumper for the filtered mesh can be directly taken from the // ElementGroup and then registered as "wheels_elements" dumper. DumperIOHelper & wheels = mesh.getGroupDumper("paraview_wheels", "wheels"); mesh.registerExternalDumper(wheels, "wheels"); mesh.setDirectoryToDumper("wheels", "./paraview/dumpable"); // Arrays and ElementTypeMapArrays can be added as external fields directly mesh.addDumpFieldExternal("displacement", displacement); ElementTypeMapArrayFilter filtered_colour( colour, wheels_elements.getElements()); - dumper::Field * colour_field_wheel = - new dumper::ElementalField(filtered_colour); + auto colour_field_wheel = + std::make_shared>( + filtered_colour); mesh.addDumpFieldExternal("color", colour_field_wheel); mesh.addDumpFieldExternalToDumper("wheels", "displacement", displacement); mesh.addDumpFieldExternalToDumper("wheels", "colour", colour); // For some specific cases the Fields should be created, as when you want to // pad an array - dumper::Field * displacement_vector_field = + auto displacement_vector_field = mesh.createNodalField(&displacement, "all", 3); mesh.addDumpFieldExternal("displacement_as_paraview_vector", displacement_vector_field); mesh.addDumpFieldExternalToDumper("wheels", "displacement_as_paraview_vector", displacement_vector_field); /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ // Fill the ElementTypeMapArray colour. fillColour(mesh, colour); /// Apply displacement and wheels rotation. Real tot_displacement = 50.; Real radius = 1.; UInt nb_steps = 100; Real theta = tot_displacement / radius; Vector l_center(spatial_dimension); Vector r_center(spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { l_center(i) = node(14, i); r_center(i) = node(2, i); } for (UInt i = 0; i < nb_steps; ++i) { displacement.clear(); Real step_ratio = Real(i) / Real(nb_steps); Real angle = step_ratio * theta; applyRotation(l_center, angle, node, displacement, lnode_1); applyRotation(l_center, angle, node, displacement, lnode_2); applyRotation(r_center, angle, node, displacement, rnode_1); applyRotation(r_center, angle, node, displacement, rnode_2); for (UInt j = 0; j < nb_nodes; ++j) { displacement(j, _x) += step_ratio * tot_displacement; } /// Dump call is finally made through Dumpable interface. mesh.dump(); mesh.dump("wheels"); } finalize(); return 0; } diff --git a/examples/io/dumper/dumper_low_level.cc b/examples/io/dumper/dumper_low_level.cc index 72d587b18..30ab30377 100644 --- a/examples/io/dumper/dumper_low_level.cc +++ b/examples/io/dumper/dumper_low_level.cc @@ -1,193 +1,194 @@ /** * @file dumper_low_level.cc * * @author Fabian Barras * @author Nicolas Richart * * @date creation: Mon Aug 17 2015 * * @brief Example of dumper::DumperIOHelper low-level methods. * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "group_manager.hh" #include "mesh.hh" #include "dumper_elemental_field.hh" #include "dumper_nodal_field.hh" #include "dumper_iohelper_paraview.hh" #include "locomotive_tools.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { /* This example aims at illustrating how to manipulate low-level methods of DumperIOHelper. The aims is to visualize a colorized moving train with Paraview */ initialize(argc, argv); // To start let us load the swiss train mesh and its mesh data information. // We aknowledge here a weel-known swiss industry for mesh donation. UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); mesh.read("swiss_train.msh"); Array & nodes = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); /* swiss_train.msh has the following physical groups that can be viewed with GMSH: "$MeshFormat 2.2 0 8 $EndMeshFormat $PhysicalNames 6 2 1 "red" 2 2 "white" 2 3 "lwheel_1" 2 4 "lwheel_2" 2 5 "rwheel_2" 2 6 "rwheel_1" $EndPhysicalNames ..." */ // Grouping nodes and elements belonging to train wheels (=four mesh data) ElementGroup & wheels_elements = mesh.createElementGroup("wheels", spatial_dimension); wheels_elements.append(mesh.getElementGroup("lwheel_1")); wheels_elements.append(mesh.getElementGroup("lwheel_2")); wheels_elements.append(mesh.getElementGroup("rwheel_1")); wheels_elements.append(mesh.getElementGroup("rwheel_2")); const Array & lnode_1 = (mesh.getElementGroup("lwheel_1")).getNodes(); const Array & lnode_2 = (mesh.getElementGroup("lwheel_2")).getNodes(); const Array & rnode_1 = (mesh.getElementGroup("rwheel_1")).getNodes(); const Array & rnode_2 = (mesh.getElementGroup("rwheel_2")).getNodes(); /* Note this Array is constructed with three components in order to warp train deformation on Paraview. A more appropriate way to do this is to set a padding in the NodalField (See example_dumpable_interface.cc.) */ Array displacement(nb_nodes, 3); // ElementalField are constructed with an ElementTypeMapArray. ElementTypeMapArray colour; colour.initialize(mesh, _with_nb_element = true); /* ------------------------------------------------------------------------ */ /* Dumper creation */ /* ------------------------------------------------------------------------ */ // Creation of two DumperParaview. One for full mesh, one for a filtered // mesh. DumperParaview dumper("train", "./paraview/dumper", false); DumperParaview wheels("wheels", "./paraview/dumper", false); // Register the full mesh dumper.registerMesh(mesh); // Register a filtered mesh limited to nodes and elements from wheels groups wheels.registerFilteredMesh(mesh, wheels_elements.getElements(), wheels_elements.getNodes()); // Generate an output file of the two mesh registered. dumper.dump(); wheels.dump(); /* At this stage no fields are attached to the two dumpers. To do so, a dumper::Field object has to be created. Several types of dumper::Field exist. In this example we present two of them. NodalField to describe nodal displacements of our train. ElementalField handling the color of our different part. */ // NodalField are constructed with an Array. - dumper::Field * displ_field = new dumper::NodalField(displacement); - dumper::Field * colour_field = new dumper::ElementalField(colour); + auto displ_field = std::make_shared>(displacement); + auto colour_field = std::make_shared>(colour); // Register the freshly created fields to our dumper. dumper.registerField("displacement", displ_field); dumper.registerField("colour", colour_field); // For the dumper wheels, fields have to be filtered at registration. // Filtered NodalField can be simply registered by adding an Array // listing the nodes. - dumper::Field * displ_field_wheel = new dumper::NodalField( + auto displ_field_wheel = std::make_shared>( displacement, 0, 0, &(wheels_elements.getNodes())); wheels.registerField("displacement", displ_field_wheel); // For the ElementalField, an ElementTypeMapArrayFilter has to be created. ElementTypeMapArrayFilter filtered_colour( colour, wheels_elements.getElements()); - dumper::Field * colour_field_wheel = - new dumper::ElementalField(filtered_colour); + auto colour_field_wheel = + std::make_shared>( + filtered_colour); wheels.registerField("colour", colour_field_wheel); /* ------------------------------------------------------------------------ */ // Now that the dumpers are created and the fields are associated, let's // paint and move the train! // Fill the ElementTypeMapArray colour according to mesh data information. fillColour(mesh, colour); // Apply displacement and wheels rotation. Real tot_displacement = 50.; Real radius = 1.; UInt nb_steps = 100; Real theta = tot_displacement / radius; Vector l_center(3); Vector r_center(3); for (UInt i = 0; i < spatial_dimension; ++i) { l_center(i) = nodes(14, i); r_center(i) = nodes(2, i); } for (UInt i = 0; i < nb_steps; ++i) { displacement.clear(); Real angle = (Real)i / (Real)nb_steps * theta; applyRotation(l_center, angle, nodes, displacement, lnode_1); applyRotation(l_center, angle, nodes, displacement, lnode_2); applyRotation(r_center, angle, nodes, displacement, rnode_1); applyRotation(r_center, angle, nodes, displacement, rnode_2); for (UInt j = 0; j < nb_nodes; ++j) { displacement(j, 0) += (Real)i / (Real)nb_steps * tot_displacement; } // Output results after each moving steps for main and wheel dumpers. dumper.dump(); wheels.dump(); } finalize(); return 0; } diff --git a/src/io/dumper/dumpable.cc b/src/io/dumper/dumpable.cc index e35791864..a840987bf 100644 --- a/src/io/dumper/dumpable.cc +++ b/src/io/dumper/dumpable.cc @@ -1,285 +1,285 @@ /** * @file dumpable.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Tue Feb 20 2018 * * @brief Implementation of the dumpable 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 "dumpable.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include namespace akantu { /* -------------------------------------------------------------------------- */ Dumpable::Dumpable() : default_dumper("") {} /* -------------------------------------------------------------------------- */ Dumpable::~Dumpable() { auto it = dumpers.begin(); auto end = dumpers.end(); for (; it != end; ++it) { delete it->second; } } /* -------------------------------------------------------------------------- */ void Dumpable::registerExternalDumper(DumperIOHelper & dumper, const std::string & dumper_name, const bool is_default) { this->dumpers[dumper_name] = &dumper; if (is_default) this->default_dumper = dumper_name; } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpMesh(const Mesh & mesh, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { this->addDumpMeshToDumper(this->default_dumper, mesh, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpMeshToDumper(const std::string & dumper_name, const Mesh & mesh, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerMesh(mesh, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpFilteredMesh( const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { this->addDumpFilteredMeshToDumper(this->default_dumper, mesh, elements_filter, nodes_filter, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpFilteredMeshToDumper( const std::string & dumper_name, const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerFilteredMesh(mesh, elements_filter, nodes_filter, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpField(const std::string & field_id) { this->addDumpFieldToDumper(this->default_dumper, field_id); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpFieldToDumper(__attribute__((unused)) const std::string & dumper_name, __attribute__((unused)) const std::string & field_id) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpFieldExternal(const std::string & field_id, - dumper::Field * field) { + std::shared_ptr field) { this->addDumpFieldExternalToDumper(this->default_dumper, field_id, field); } /* -------------------------------------------------------------------------- */ -void Dumpable::addDumpFieldExternalToDumper(const std::string & dumper_name, - const std::string & field_id, - dumper::Field * field) { +void Dumpable::addDumpFieldExternalToDumper( + const std::string & dumper_name, const std::string & field_id, + std::shared_ptr field) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerField(field_id, field); } /* -------------------------------------------------------------------------- */ void Dumpable::removeDumpField(const std::string & field_id) { this->removeDumpFieldFromDumper(this->default_dumper, field_id); } /* -------------------------------------------------------------------------- */ void Dumpable::removeDumpFieldFromDumper(const std::string & dumper_name, const std::string & field_id) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.unRegisterField(field_id); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpFieldVector(const std::string & field_id) { this->addDumpFieldVectorToDumper(this->default_dumper, field_id); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpFieldVectorToDumper(__attribute__((unused)) const std::string & dumper_name, __attribute__((unused)) const std::string & field_id) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpFieldTensor(const std::string & field_id) { this->addDumpFieldTensorToDumper(this->default_dumper, field_id); } /* -------------------------------------------------------------------------- */ void Dumpable::addDumpFieldTensorToDumper(__attribute__((unused)) const std::string & dumper_name, __attribute__((unused)) const std::string & field_id) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void Dumpable::setDirectory(const std::string & directory) { this->setDirectoryToDumper(this->default_dumper, directory); } /* -------------------------------------------------------------------------- */ void Dumpable::setDirectoryToDumper(const std::string & dumper_name, const std::string & directory) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Dumpable::setBaseName(const std::string & basename) { this->setBaseNameToDumper(this->default_dumper, basename); } /* -------------------------------------------------------------------------- */ void Dumpable::setBaseNameToDumper(const std::string & dumper_name, const std::string & basename) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.setBaseName(basename); } /* -------------------------------------------------------------------------- */ void Dumpable::setTimeStepToDumper(Real time_step) { this->setTimeStepToDumper(this->default_dumper, time_step); } /* -------------------------------------------------------------------------- */ void Dumpable::setTimeStepToDumper(const std::string & dumper_name, Real time_step) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void Dumpable::setTextModeToDumper(const std::string & dumper_name) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.getDumper().setMode(iohelper::TEXT); } /* -------------------------------------------------------------------------- */ void Dumpable::setTextModeToDumper() { DumperIOHelper & dumper = this->getDumper(this->default_dumper); dumper.getDumper().setMode(iohelper::TEXT); } /* -------------------------------------------------------------------------- */ void Dumpable::dump(const std::string & dumper_name) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.dump(); } /* -------------------------------------------------------------------------- */ void Dumpable::dump() { this->dump(this->default_dumper); } /* -------------------------------------------------------------------------- */ void Dumpable::dump(const std::string & dumper_name, UInt step) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.dump(step); } /* -------------------------------------------------------------------------- */ void Dumpable::dump(UInt step) { this->dump(this->default_dumper, step); } /* -------------------------------------------------------------------------- */ void Dumpable::dump(const std::string & dumper_name, Real time, UInt step) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.dump(time, step); } /* -------------------------------------------------------------------------- */ void Dumpable::dump(Real time, UInt step) { this->dump(this->default_dumper, time, step); } /* -------------------------------------------------------------------------- */ -void Dumpable::internalAddDumpFieldToDumper(const std::string & dumper_name, - const std::string & field_id, - dumper::Field * field) { +void Dumpable::internalAddDumpFieldToDumper( + const std::string & dumper_name, const std::string & field_id, + std::shared_ptr field) { DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerField(field_id, field); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Dumpable::getDumper() { return this->getDumper(this->default_dumper); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Dumpable::getDumper(const std::string & dumper_name) { auto it = this->dumpers.find(dumper_name); auto end = this->dumpers.end(); if (it == end) AKANTU_EXCEPTION("Dumper " << dumper_name << "has not been registered, yet."); return *(it->second); } /* -------------------------------------------------------------------------- */ std::string Dumpable::getDefaultDumperName() const { return this->default_dumper; } -} // akantu +} // namespace akantu #endif diff --git a/src/io/dumper/dumpable_dummy.hh b/src/io/dumper/dumpable_dummy.hh index eb9420dfd..f60256253 100644 --- a/src/io/dumper/dumpable_dummy.hh +++ b/src/io/dumper/dumpable_dummy.hh @@ -1,264 +1,265 @@ /** * @file dumpable_dummy.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Oct 26 2012 * @date last modification: Tue Feb 20 2018 * * @brief Interface for object who wants to dump themselves * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DUMPABLE_DUMMY_HH__ #define __AKANTU_DUMPABLE_DUMMY_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused" namespace dumper { class Field; } class DumperIOHelper; class Mesh; /* -------------------------------------------------------------------------- */ class Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Dumpable(){}; virtual ~Dumpable(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template inline void registerDumper(const std::string & dumper_name, const std::string & file_name = "", const bool is_default = false) {} void registerExternalDumper(DumperIOHelper * dumper, const std::string & dumper_name, const bool is_default = false) {} void addDumpMesh(const Mesh & mesh, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined) {} void addDumpMeshToDumper(const std::string & dumper_name, const Mesh & mesh, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined) { } void addDumpFilteredMesh(const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined) { } void addDumpFilteredMeshToDumper( const std::string & dumper_name, const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined) {} virtual void addDumpField(const std::string & field_id) { AKANTU_TO_IMPLEMENT(); } virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_TO_IMPLEMENT(); } virtual void addDumpFieldExternal(const std::string & field_id, - dumper::Field * field) { + std::shared_ptr field) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } - virtual void addDumpFieldExternalToDumper(const std::string & dumper_name, - const std::string & field_id, - dumper::Field * field) { + virtual void + addDumpFieldExternalToDumper(const std::string & dumper_name, + const std::string & field_id, + std::shared_ptr field) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } template void addDumpFieldExternal(const std::string & field_id, const Array & field) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } template void addDumpFieldExternalToDumper(const std::string & dumper_name, const std::string & field_id, const Array & field) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } template void addDumpFieldExternal(const std::string & field_id, const ElementTypeMapArray & field, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } template void addDumpFieldExternalToDumper( const std::string & dumper_name, const std::string & field_id, const ElementTypeMapArray & field, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void removeDumpField(const std::string & field_id) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void removeDumpFieldFromDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void setDirecory(const std::string & directory) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void setDirectoryToDumper(const std::string & dumper_name, const std::string & directory) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void setBaseName(const std::string & basename) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void setBaseNameToDumper(const std::string & dumper_name, const std::string & basename) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void setTextModeToDumper(const std::string & dumper_name) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void setTextModeToDumper() { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void dump() { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void dump(const std::string & dumper_name) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void dump(UInt step) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void dump(const std::string & dumper_name, UInt step) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void dump(Real current_time, UInt step) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } void dump(const std::string & dumper_name, Real current_time, UInt step) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } protected: void internalAddDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id, - dumper::Field * field) { + std::shared_ptr field) { AKANTU_DEBUG_WARNING("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } protected: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: DumperIOHelper & getDumper() { AKANTU_ERROR("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } DumperIOHelper & getDumper(const std::string & dumper_name) { AKANTU_ERROR("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } template T & getDumper(const std::string & dumper_name) { AKANTU_ERROR("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } std::string getDefaultDumperName() { AKANTU_ERROR("No dumper activated at compilation, turn on " "AKANTU_USE_IOHELPER in cmake."); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: }; #pragma GCC diagnostic pop -} // akantu +} // namespace akantu #endif /* __AKANTU_DUMPABLE_DUMMY_HH__ */ diff --git a/src/io/dumper/dumpable_inline_impl.hh b/src/io/dumper/dumpable_inline_impl.hh index 012b570b4..822b5b957 100644 --- a/src/io/dumper/dumpable_inline_impl.hh +++ b/src/io/dumper/dumpable_inline_impl.hh @@ -1,134 +1,134 @@ /** * @file dumpable_inline_impl.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Nov 08 2017 * * @brief Implementation of the Dumpable class * * @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_DUMPABLE_INLINE_IMPL_HH__ #define __AKANTU_DUMPABLE_INLINE_IMPL_HH__ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumper_elemental_field.hh" #include "dumper_nodal_field.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template inline void Dumpable::registerDumper(const std::string & dumper_name, const std::string & file_name, const bool is_default) { if (this->dumpers.find(dumper_name) != this->dumpers.end()) { AKANTU_DEBUG_INFO("Dumper " + dumper_name + "is already registered."); } std::string name = file_name; if (name == "") name = dumper_name; this->dumpers[dumper_name] = new T(name); if (is_default) this->default_dumper = dumper_name; } /* -------------------------------------------------------------------------- */ template inline void Dumpable::addDumpFieldExternal(const std::string & field_id, const Array & field) { this->addDumpFieldExternalToDumper(this->default_dumper, field_id, field); } /* -------------------------------------------------------------------------- */ template inline void Dumpable::addDumpFieldExternalToDumper(const std::string & dumper_name, const std::string & field_id, const Array & field) { - dumper::Field * field_cont = new dumper::NodalField(field); + auto field_cont = std::make_shared>(field); DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerField(field_id, field_cont); } /* -------------------------------------------------------------------------- */ template inline void Dumpable::addDumpFieldExternal(const std::string & field_id, const ElementTypeMapArray & field, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { this->addDumpFieldExternalToDumper(this->default_dumper, field_id, field, spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ template inline void Dumpable::addDumpFieldExternalToDumper( const std::string & dumper_name, const std::string & field_id, const ElementTypeMapArray & field, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { - dumper::Field * field_cont; + std::shared_ptr field_cont; #if defined(AKANTU_IGFEM) if (element_kind == _ek_igfem) { - field_cont = new dumper::IGFEMElementalField(field, spatial_dimension, - ghost_type, element_kind); + field_cont = std::make_shared>( + field, spatial_dimension, ghost_type, element_kind); } else #endif - field_cont = new dumper::ElementalField(field, spatial_dimension, - ghost_type, element_kind); + field_cont = std::make_shared>( + field, spatial_dimension, ghost_type, element_kind); DumperIOHelper & dumper = this->getDumper(dumper_name); dumper.registerField(field_id, field_cont); } /* -------------------------------------------------------------------------- */ template inline T & Dumpable::getDumper(const std::string & dumper_name) { DumperIOHelper & dumper = this->getDumper(dumper_name); try { auto & templated_dumper = dynamic_cast(dumper); return templated_dumper; } catch (std::bad_cast &) { AKANTU_EXCEPTION("Dumper " << dumper_name << " is not of type: " << debug::demangle(typeid(T).name())); } } /* -------------------------------------------------------------------------- */ -} // akantu +} // namespace akantu #endif #endif /* __AKANTU_DUMPABLE_INLINE_IMPL_HH__ */ diff --git a/src/io/dumper/dumpable_iohelper.hh b/src/io/dumper/dumpable_iohelper.hh index 02bacd64f..5a5995951 100644 --- a/src/io/dumper/dumpable_iohelper.hh +++ b/src/io/dumper/dumpable_iohelper.hh @@ -1,192 +1,193 @@ /** * @file dumpable_iohelper.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Tue Jan 06 2015 * @date last modification: Sun Dec 03 2017 * * @brief Interface for object who wants to dump themselves * * @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 "dumper_iohelper.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DUMPABLE_IOHELPER_HH__ #define __AKANTU_DUMPABLE_IOHELPER_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Dumpable(); virtual ~Dumpable(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// create a new dumper (of templated type T) and register it under /// dumper_name. file_name is used for construction of T. is default states if /// this dumper is the default dumper. template inline void registerDumper(const std::string & dumper_name, const std::string & file_name = "", const bool is_default = false); /// register an externally created dumper void registerExternalDumper(DumperIOHelper & dumper, const std::string & dumper_name, const bool is_default = false); /// register a mesh to the default dumper void addDumpMesh(const Mesh & mesh, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); /// register a mesh to the default identified by its name void addDumpMeshToDumper(const std::string & dumper_name, const Mesh & mesh, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); /// register a filtered mesh as the default dumper void addDumpFilteredMesh(const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); /// register a filtered mesh and provides a name void addDumpFilteredMeshToDumper( const std::string & dumper_name, const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); /// to implement virtual void addDumpField(const std::string & field_id); /// to implement virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /// add a field virtual void addDumpFieldExternal(const std::string & field_id, - dumper::Field * field); - virtual void addDumpFieldExternalToDumper(const std::string & dumper_name, - const std::string & field_id, - dumper::Field * field); + std::shared_ptr field); + virtual void + addDumpFieldExternalToDumper(const std::string & dumper_name, + const std::string & field_id, + std::shared_ptr field); template inline void addDumpFieldExternal(const std::string & field_id, const Array & field); template inline void addDumpFieldExternalToDumper(const std::string & dumper_name, const std::string & field_id, const Array & field); template inline void addDumpFieldExternal(const std::string & field_id, const ElementTypeMapArray & field, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); template inline void addDumpFieldExternalToDumper( const std::string & dumper_name, const std::string & field_id, const ElementTypeMapArray & field, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); void removeDumpField(const std::string & field_id); void removeDumpFieldFromDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldVector(const std::string & field_id); virtual void addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensor(const std::string & field_id); virtual void addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id); void setDirectory(const std::string & directory); void setDirectoryToDumper(const std::string & dumper_name, const std::string & directory); void setBaseName(const std::string & basename); void setBaseNameToDumper(const std::string & dumper_name, const std::string & basename); void setTimeStepToDumper(Real time_step); void setTimeStepToDumper(const std::string & dumper_name, Real time_step); void setTextModeToDumper(const std::string & dumper_name); void setTextModeToDumper(); virtual void dump(); virtual void dump(UInt step); virtual void dump(Real time, UInt step); virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); public: void internalAddDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id, - dumper::Field * field); + std::shared_ptr field); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: DumperIOHelper & getDumper(); DumperIOHelper & getDumper(const std::string & dumper_name); template T & getDumper(const std::string & dumper_name); std::string getDefaultDumperName() const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: using DumperMap = std::map; using DumperSet = std::set; DumperMap dumpers; std::string default_dumper; }; -} // akantu +} // namespace akantu #endif /* __AKANTU_DUMPABLE_IOHELPER_HH__ */ diff --git a/src/io/dumper/dumper_compute.hh b/src/io/dumper/dumper_compute.hh index f16c1c9e5..ff8fa5a46 100644 --- a/src/io/dumper/dumper_compute.hh +++ b/src/io/dumper/dumper_compute.hh @@ -1,268 +1,270 @@ /** * @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 auto & old_nb_components = this->sub_field.getNbComponents(dim, ghost_type, kind); - for(auto type : old_nb_components.elementTypes(dim, ghost_type, kind)) { + 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; + inline std::shared_ptr connect(FieldComputeProxy & proxy) override; /// for connection to a FieldCompute - ComputeFunctorInterface * connect(HomogenizerProxy & proxy) override; + std::shared_ptr + 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) { + inline static std::shared_ptr + createFieldCompute(std::shared_ptr 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) { + template std::shared_ptr 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 + std::shared_ptr connectToFunctor(T * ptr) { + return std::make_shared>(*ptr, func); } template - Field * + std::shared_ptr connectToFunctor(__attribute__((unused)) FieldCompute, return_type2> * ptr) { throw; // return new FieldCompute(*ptr,func); return nullptr; } template - Field * connectToFunctor( + std::shared_ptr 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 * +inline std::shared_ptr FieldCompute::connect(FieldComputeProxy & proxy) { return proxy.connectToField(this); } /* -------------------------------------------------------------------------- */ __END_AKANTU_DUMPER__ -} // akantu +} // namespace akantu #endif /* __AKANTU_DUMPER_COMPUTE_HH__ */ diff --git a/src/io/dumper/dumper_field.hh b/src/io/dumper/dumper_field.hh index ebe765ce8..5a2c58491 100644 --- a/src/io/dumper/dumper_field.hh +++ b/src/io/dumper/dumper_field.hh @@ -1,137 +1,137 @@ /** * @file dumper_field.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Sep 02 2014 * @date last modification: Tue Feb 20 2018 * * @brief Common interface for 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_FIELD_HH__ #define __AKANTU_DUMPER_FIELD_HH__ /* -------------------------------------------------------------------------- */ #include "dumper_iohelper.hh" /* -------------------------------------------------------------------------- */ namespace akantu { __BEGIN_AKANTU_DUMPER__ /* -------------------------------------------------------------------------- */ class FieldComputeProxy; class FieldComputeBaseInterface; class ComputeFunctorInterface; class HomogenizerProxy; /* -------------------------------------------------------------------------- */ /// Field interface class Field { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Field() = default; virtual ~Field() = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: #ifdef AKANTU_USE_IOHELPER /// register this to the provided dumper virtual void registerToDumper(const std::string & id, iohelper::Dumper & dumper) = 0; #endif /// set the number of data per item (used for elements fields at the moment) virtual void setNbData(__attribute__((unused)) UInt nb_data) { AKANTU_TO_IMPLEMENT(); }; /// set the number of data per elem (used for elements fields at the moment) virtual void setNbDataPerElem(__attribute__((unused)) const ElementTypeMap & nb_data) { AKANTU_TO_IMPLEMENT(); }; /// set the number of data per elem (used for elements fields at the moment) virtual void setNbDataPerElem(__attribute__((unused)) UInt nb_data) { AKANTU_TO_IMPLEMENT(); }; /// get the number of components of the hosted field virtual ElementTypeMap getNbComponents(__attribute__((unused)) UInt dim = _all_dimensions, __attribute__((unused)) GhostType ghost_type = _not_ghost, __attribute__((unused)) ElementKind kind = _ek_not_defined) { throw; }; /// for connection to a FieldCompute - inline virtual Field * connect(__attribute__((unused)) - FieldComputeProxy & proxy) { + inline virtual std::shared_ptr connect(__attribute__((unused)) + FieldComputeProxy & proxy) { throw; }; /// for connection to a FieldCompute - inline virtual ComputeFunctorInterface * connect(__attribute__((unused)) - HomogenizerProxy & proxy) { + inline virtual std::shared_ptr + connect(__attribute__((unused)) HomogenizerProxy & proxy) { throw; }; /// check if the same quantity of data for all element types virtual void checkHomogeneity() = 0; /// return the dumper name std::string getGroupName() { return group_name; }; /// return the id of the field std::string getID() { return field_id; }; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the flag to know if the field is homogeneous/contiguous virtual bool isHomogeneous() { return homogeneous; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the flag to know if it is homogeneous bool homogeneous{false}; /// the name of the group it was associated to std::string group_name; /// the name of the dumper it was associated to std::string field_id; }; /* -------------------------------------------------------------------------- */ __END_AKANTU_DUMPER__ -} // akantu +} // namespace akantu #endif /* __AKANTU_DUMPER_FIELD_HH__ */ diff --git a/src/io/dumper/dumper_generic_elemental_field.hh b/src/io/dumper/dumper_generic_elemental_field.hh index a66095318..4d4c5186a 100644 --- a/src/io/dumper/dumper_generic_elemental_field.hh +++ b/src/io/dumper/dumper_generic_elemental_field.hh @@ -1,215 +1,218 @@ /** * @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 { + inline std::shared_ptr connect(FieldComputeProxy & proxy) override { return proxy.connectToField(this); } /// for connection to a Homogenizer - inline ComputeFunctorInterface * connect(HomogenizerProxy & proxy) override { + inline std::shared_ptr + connect(HomogenizerProxy & proxy) override { return proxy.connectToField(this); }; virtual iterator begin() { /// type iterators on the elemental field - auto types = this->field.elementTypes(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 for (; tit != end && this->field(*tit, this->ghost_type).size() == 0; ++tit) { } auto type = *tit; if (tit == end) return this->end(); /// getting information for the field of the given type const auto & vect = this->field(type, this->ghost_type); UInt nb_data_per_elem = this->getNbDataPerElem(type); /// define element-wise iterator 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() { - auto types = this->field.elementTypes(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(); 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); 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) { - auto tit = this->field.elementTypes( - this->spatial_dimension, this->ghost_type, this->element_kind).begin(); + 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 +} // namespace akantu #endif /* __AKANTU_DUMPER_GENERIC_ELEMENTAL_FIELD_HH__ */ diff --git a/src/io/dumper/dumper_homogenizing_field.hh b/src/io/dumper/dumper_homogenizing_field.hh index daf4c9eb9..ae12276cc 100644 --- a/src/io/dumper/dumper_homogenizing_field.hh +++ b/src/io/dumper/dumper_homogenizing_field.hh @@ -1,216 +1,215 @@ /** * @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(const Vector & , - Vector & res, +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) { 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); + inline static std::shared_ptr + createHomogenizer(Field & field); template - inline ComputeFunctorInterface * connectToField(T * field) { + inline std::shared_ptr connectToField(T * field) { ElementTypeMap nb_components = field->getNbComponents(); using ret_type = typename T::types::return_type; return this->instantiateHomogenizer(nb_components); } template - inline ComputeFunctorInterface * + inline std::shared_ptr instantiateHomogenizer(ElementTypeMap & nb_components); }; /* -------------------------------------------------------------------------- */ template -inline ComputeFunctorInterface * +inline std::shared_ptr HomogenizerProxy::instantiateHomogenizer(ElementTypeMap & nb_components) { using Homogenizer = dumper::AvgHomogenizingFunctor; - auto * foo = new Homogenizer(nb_components); - return foo; + return std::make_shared(nb_components); } template <> -inline ComputeFunctorInterface * +inline std::shared_ptr HomogenizerProxy::instantiateHomogenizer>( __attribute__((unused)) ElementTypeMap & nb_components) { throw; return nullptr; } /* -------------------------------------------------------------------------- */ /// for connection to a FieldCompute template -inline ComputeFunctorInterface * +inline std::shared_ptr FieldCompute::connect(HomogenizerProxy & proxy) { return proxy.connectToField(this); } /* -------------------------------------------------------------------------- */ -inline ComputeFunctorInterface * +inline std::shared_ptr 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/dumper/dumper_iohelper.cc b/src/io/dumper/dumper_iohelper.cc index c963934e1..7ca4dec82 100644 --- a/src/io/dumper/dumper_iohelper.cc +++ b/src/io/dumper/dumper_iohelper.cc @@ -1,319 +1,314 @@ /** * @file dumper_iohelper.cc * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Oct 26 2012 * @date last modification: Tue Feb 20 2018 * * @brief implementation of DumperIOHelper * * @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 "dumper_elemental_field.hh" #include "dumper_filtered_connectivity.hh" #include "dumper_iohelper.hh" #include "dumper_nodal_field.hh" #include "dumper_variable.hh" #include "mesh.hh" #if defined(AKANTU_IGFEM) #include "dumper_igfem_connectivity.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DumperIOHelper::DumperIOHelper() = default; /* -------------------------------------------------------------------------- */ -DumperIOHelper::~DumperIOHelper() { - for (auto it = fields.begin(); it != fields.end(); ++it) { - delete it->second; - } - - delete dumper; -} +DumperIOHelper::~DumperIOHelper() {} /* -------------------------------------------------------------------------- */ void DumperIOHelper::setParallelContext(bool is_parallel) { UInt whoami = Communicator::getStaticCommunicator().whoAmI(); UInt nb_proc = Communicator::getStaticCommunicator().getNbProc(); if (is_parallel) dumper->setParallelContext(whoami, nb_proc); else dumper->setParallelContext(0, 1); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::setDirectory(const std::string & directory) { this->directory = directory; dumper->setPrefix(directory); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::setBaseName(const std::string & basename) { filename = basename; } /* -------------------------------------------------------------------------- */ void DumperIOHelper::setTimeStep(Real time_step) { if (!time_activated) this->dumper->activateTimeDescFiles(time_step); else this->dumper->setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::dump() { try { dumper->dump(filename, count); } catch (iohelper::IOHelperException & e) { AKANTU_ERROR( "I was not able to dump your data with a Dumper: " << e.what()); } ++count; } /* -------------------------------------------------------------------------- */ void DumperIOHelper::dump(UInt step) { this->count = step; this->dump(); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::dump(Real current_time, UInt step) { this->dumper->setCurrentTime(current_time); this->dump(step); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::registerMesh(const Mesh & mesh, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { #if defined(AKANTU_IGFEM) if (element_kind == _ek_igfem) { registerField("connectivities", new dumper::IGFEMConnectivityField( mesh.getConnectivities(), spatial_dimension, ghost_type)); } else #endif registerField("connectivities", - new dumper::ElementalField(mesh.getConnectivities(), - spatial_dimension, - ghost_type, element_kind)); + std::make_shared>( + mesh.getConnectivities(), spatial_dimension, ghost_type, + element_kind)); - registerField("positions", new dumper::NodalField(mesh.getNodes())); + registerField("positions", + std::make_shared>(mesh.getNodes())); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::registerFilteredMesh( const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) { auto * f_connectivities = new ElementTypeMapArrayFilter( mesh.getConnectivities(), elements_filter); this->registerField("connectivities", - new dumper::FilteredConnectivityField( + std::make_shared( *f_connectivities, nodes_filter, spatial_dimension, ghost_type, element_kind)); - this->registerField("positions", new dumper::NodalField( - mesh.getNodes(), 0, 0, &nodes_filter)); + this->registerField("positions", + std::make_shared>( + mesh.getNodes(), 0, 0, &nodes_filter)); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::registerField(const std::string & field_id, - dumper::Field * field) { + std::shared_ptr field) { auto it = fields.find(field_id); if (it != fields.end()) { AKANTU_DEBUG_WARNING( "The field " << field_id << " is already registered in this Dumper. Field ignored."); return; } fields[field_id] = field; field->registerToDumper(field_id, *dumper); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::unRegisterField(const std::string & field_id) { auto it = fields.find(field_id); if (it == fields.end()) { AKANTU_DEBUG_WARNING( "The field " << field_id << " is not registered in this Dumper. Nothing to do."); return; } - delete it->second; fields.erase(it); } /* -------------------------------------------------------------------------- */ -void DumperIOHelper::registerVariable(const std::string & variable_id, - dumper::VariableBase * variable) { +void DumperIOHelper::registerVariable( + const std::string & variable_id, + std::shared_ptr variable) { auto it = variables.find(variable_id); if (it != variables.end()) { AKANTU_DEBUG_WARNING( "The Variable " << variable_id << " is already registered in this Dumper. Variable ignored."); return; } variables[variable_id] = variable; variable->registerToDumper(variable_id, *dumper); } /* -------------------------------------------------------------------------- */ void DumperIOHelper::unRegisterVariable(const std::string & variable_id) { auto it = variables.find(variable_id); if (it == variables.end()) { AKANTU_DEBUG_WARNING( "The variable " << variable_id << " is not registered in this Dumper. Nothing to do."); return; } - delete it->second; variables.erase(it); } /* -------------------------------------------------------------------------- */ template iohelper::ElemType getIOHelperType() { AKANTU_TO_IMPLEMENT(); return iohelper::MAX_ELEM_TYPE; } template <> iohelper::ElemType getIOHelperType<_point_1>() { return iohelper::POINT_SET; } template <> iohelper::ElemType getIOHelperType<_segment_2>() { return iohelper::LINE1; } template <> iohelper::ElemType getIOHelperType<_segment_3>() { return iohelper::LINE2; } template <> iohelper::ElemType getIOHelperType<_triangle_3>() { return iohelper::TRIANGLE1; } template <> iohelper::ElemType getIOHelperType<_triangle_6>() { return iohelper::TRIANGLE2; } template <> iohelper::ElemType getIOHelperType<_quadrangle_4>() { return iohelper::QUAD1; } template <> iohelper::ElemType getIOHelperType<_quadrangle_8>() { return iohelper::QUAD2; } template <> iohelper::ElemType getIOHelperType<_tetrahedron_4>() { return iohelper::TETRA1; } template <> iohelper::ElemType getIOHelperType<_tetrahedron_10>() { return iohelper::TETRA2; } template <> iohelper::ElemType getIOHelperType<_hexahedron_8>() { return iohelper::HEX1; } template <> iohelper::ElemType getIOHelperType<_hexahedron_20>() { return iohelper::HEX2; } template <> iohelper::ElemType getIOHelperType<_pentahedron_6>() { return iohelper::PRISM1; } template <> iohelper::ElemType getIOHelperType<_pentahedron_15>() { return iohelper::PRISM2; } #if defined(AKANTU_COHESIVE_ELEMENT) template <> iohelper::ElemType getIOHelperType<_cohesive_1d_2>() { return iohelper::COH1D2; } template <> iohelper::ElemType getIOHelperType<_cohesive_2d_4>() { return iohelper::COH2D4; } template <> iohelper::ElemType getIOHelperType<_cohesive_2d_6>() { return iohelper::COH2D6; } template <> iohelper::ElemType getIOHelperType<_cohesive_3d_6>() { return iohelper::COH3D6; } template <> iohelper::ElemType getIOHelperType<_cohesive_3d_12>() { return iohelper::COH3D12; } template <> iohelper::ElemType getIOHelperType<_cohesive_3d_8>() { return iohelper::COH3D8; } // template <> // iohelper::ElemType getIOHelperType<_cohesive_3d_16>() { return // iohelper::COH3D16; } #endif #if defined(AKANTU_STRUCTURAL_MECHANICS) template <> iohelper::ElemType getIOHelperType<_bernoulli_beam_2>() { return iohelper::BEAM2; } template <> iohelper::ElemType getIOHelperType<_bernoulli_beam_3>() { return iohelper::BEAM3; } #endif /* -------------------------------------------------------------------------- */ UInt getIOHelperType(ElementType type) { UInt ioh_type = iohelper::MAX_ELEM_TYPE; #define GET_IOHELPER_TYPE(type) ioh_type = getIOHelperType(); AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_IOHELPER_TYPE); #undef GET_IOHELPER_TYPE return ioh_type; } /* -------------------------------------------------------------------------- */ -} // akantu +} // namespace akantu namespace iohelper { template <> DataType getDataType() { return getDataType>(); } -} +} // namespace iohelper diff --git a/src/io/dumper/dumper_iohelper.hh b/src/io/dumper/dumper_iohelper.hh index f824215bd..6125a253e 100644 --- a/src/io/dumper/dumper_iohelper.hh +++ b/src/io/dumper/dumper_iohelper.hh @@ -1,158 +1,160 @@ /** * @file dumper_iohelper.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Oct 26 2012 * @date last modification: Sun Dec 03 2017 * * @brief Define the akantu dumper interface for IOhelper dumpers * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_common.hh" #include "aka_types.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DUMPER_IOHELPER_HH__ #define __AKANTU_DUMPER_IOHELPER_HH__ /* -------------------------------------------------------------------------- */ namespace iohelper { class Dumper; } namespace akantu { UInt getIOHelperType(ElementType type); namespace dumper { class Field; class VariableBase; -} +} // namespace dumper class Mesh; class DumperIOHelper { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DumperIOHelper(); virtual ~DumperIOHelper(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register a given Mesh for the current dumper virtual void registerMesh(const Mesh & mesh, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); /// register a filtered Mesh (provided filter lists) for the current dumper virtual void registerFilteredMesh(const Mesh & mesh, const ElementTypeMapArray & elements_filter, const Array & nodes_filter, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_not_defined); /// register a Field object identified by name and provided by pointer - void registerField(const std::string & field_id, dumper::Field * field); + void registerField(const std::string & field_id, + std::shared_ptr field); /// remove the Field identified by name from managed fields void unRegisterField(const std::string & field_id); /// register a VariableBase object identified by name and provided by pointer void registerVariable(const std::string & variable_id, - dumper::VariableBase * variable); + std::shared_ptr variable); /// remove a VariableBase identified by name from managed fields void unRegisterVariable(const std::string & variable_id); /// request dump: this calls IOHelper dump routine virtual void dump(); /// request dump: this first set the current step and then calls IOHelper dump /// routine virtual void dump(UInt step); /// request dump: this first set the current step and current time and then /// calls IOHelper dump routine virtual void dump(Real current_time, UInt step); /// set the parallel context for IOHeper virtual void setParallelContext(bool is_parallel); /// set the directory where to generate the dumped files virtual void setDirectory(const std::string & directory); /// set the base name (needed by most IOHelper dumpers) virtual void setBaseName(const std::string & basename); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// direct access to the iohelper::Dumper object AKANTU_GET_MACRO(Dumper, *dumper, iohelper::Dumper &) /// set the timestep of the iohelper::Dumper void setTimeStep(Real time_step); public: /* ------------------------------------------------------------------------ */ /* Variable wrapper */ template ::value> class Variable; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// internal iohelper::Dumper - iohelper::Dumper * dumper; + std::unique_ptr dumper; - using Fields = std::map; - using Variables = std::map; + using Fields = std::map>; + using Variables = + std::map>; /// list of registered fields to dump Fields fields; Variables variables; /// dump counter UInt count{0}; /// directory name std::string directory; /// filename prefix std::string filename; /// is time tracking activated in the dumper bool time_activated{false}; }; -} // akantu +} // namespace akantu #endif /* __AKANTU_DUMPER_IOHELPER_HH__ */ diff --git a/src/io/dumper/dumper_iohelper_paraview.cc b/src/io/dumper/dumper_iohelper_paraview.cc index fb8531321..7a7525b94 100644 --- a/src/io/dumper/dumper_iohelper_paraview.cc +++ b/src/io/dumper/dumper_iohelper_paraview.cc @@ -1,66 +1,65 @@ /** * @file dumper_iohelper_paraview.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Sun Sep 26 2010 * @date last modification: Mon Jan 22 2018 * * @brief implementations of DumperParaview * * @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 "dumper_iohelper_paraview.hh" #include "communicator.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { DumperParaview::DumperParaview(const std::string & filename, const std::string & directory, bool parallel) : DumperIOHelper() { - iohelper::DumperParaview * dumper_para = new iohelper::DumperParaview(); - dumper = dumper_para; + dumper = std::make_unique(); setBaseName(filename); this->setParallelContext(parallel); - dumper_para->setMode(iohelper::BASE64); - dumper_para->setPrefix(directory); - dumper_para->init(); + dumper->setMode(iohelper::BASE64); + dumper->setPrefix(directory); + dumper->init(); } /* -------------------------------------------------------------------------- */ DumperParaview::~DumperParaview() = default; /* -------------------------------------------------------------------------- */ void DumperParaview::setBaseName(const std::string & basename) { DumperIOHelper::setBaseName(basename); - static_cast(dumper)->setVTUSubDirectory(filename + - "-VTU"); + static_cast(dumper.get()) + ->setVTUSubDirectory(filename + "-VTU"); } } // namespace akantu diff --git a/src/io/dumper/dumper_text.cc b/src/io/dumper/dumper_text.cc index e379a6683..cdac9052c 100644 --- a/src/io/dumper/dumper_text.cc +++ b/src/io/dumper/dumper_text.cc @@ -1,112 +1,113 @@ /** * @file dumper_text.cc * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Nov 07 2017 * * @brief implementation of text dumper * * @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 "dumper_text.hh" #include "communicator.hh" #include "dumper_nodal_field.hh" #include "mesh.hh" #include namespace akantu { /* -------------------------------------------------------------------------- */ DumperText::DumperText(const std::string & basename, iohelper::TextDumpMode mode, bool parallel) : DumperIOHelper() { AKANTU_DEBUG_IN(); - iohelper::DumperText * dumper_text = new iohelper::DumperText(mode); - this->dumper = dumper_text; + this->dumper = std::make_unique(mode); this->setBaseName(basename); this->setParallelContext(parallel); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DumperText::registerMesh(const Mesh & mesh, __attribute__((unused)) UInt spatial_dimension, __attribute__((unused)) const GhostType & ghost_type, __attribute__((unused)) const ElementKind & element_kind) { - registerField("position", new dumper::NodalField(mesh.getNodes())); + registerField("position", + std::make_shared>(mesh.getNodes())); // in parallel we need node type UInt nb_proc = mesh.getCommunicator().getNbProc(); if (nb_proc > 1) { - registerField("nodes_type", - new dumper::NodalField(mesh.getNodesFlags())); + registerField("nodes_type", std::make_shared>( + mesh.getNodesFlags())); } } /* -------------------------------------------------------------------------- */ void DumperText::registerFilteredMesh( const Mesh & mesh, __attribute__((unused)) const ElementTypeMapArray & elements_filter, const Array & nodes_filter, __attribute__((unused)) UInt spatial_dimension, __attribute__((unused)) const GhostType & ghost_type, __attribute__((unused)) const ElementKind & element_kind) { - registerField("position", new dumper::NodalField( + registerField("position", std::make_shared>( mesh.getNodes(), 0, 0, &nodes_filter)); // in parallel we need node type UInt nb_proc = mesh.getCommunicator().getNbProc(); if (nb_proc > 1) { - registerField("nodes_type", new dumper::NodalField( - mesh.getNodesFlags(), 0, 0, &nodes_filter)); + registerField("nodes_type", + std::make_shared>( + mesh.getNodesFlags(), 0, 0, &nodes_filter)); } } /* -------------------------------------------------------------------------- */ void DumperText::setBaseName(const std::string & basename) { AKANTU_DEBUG_IN(); DumperIOHelper::setBaseName(basename); - static_cast(this->dumper) + static_cast(this->dumper.get()) ->setDataSubDirectory(this->filename + "-DataFiles"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DumperText::setPrecision(UInt prec) { AKANTU_DEBUG_IN(); - static_cast(this->dumper)->setPrecision(prec); + static_cast(this->dumper.get())->setPrecision(prec); AKANTU_DEBUG_OUT(); } -} // akantu +} // namespace akantu diff --git a/src/mesh/group_manager.hh b/src/mesh/group_manager.hh index d255c3761..7c5910d58 100644 --- a/src/mesh/group_manager.hh +++ b/src/mesh/group_manager.hh @@ -1,312 +1,313 @@ /** * @file group_manager.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Nov 13 2013 * @date last modification: Wed Feb 07 2018 * * @brief Stores information relevent to the notion of element and nodes * groups. * * @section LICENSE * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GROUP_MANAGER_HH__ #define __AKANTU_GROUP_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_iterators.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class ElementGroup; class NodeGroup; class Mesh; class Element; class ElementSynchronizer; template class CommunicationBufferTemplated; namespace dumper { class Field; } } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class GroupManager { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: using ElementGroups = std::map; using NodeGroups = std::map; public: using GroupManagerTypeSet = std::set; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GroupManager(const Mesh & mesh, const ID & id = "group_manager", const MemoryID & memory_id = 0); virtual ~GroupManager(); /* ------------------------------------------------------------------------ */ /* Groups iterators */ /* ------------------------------------------------------------------------ */ public: using node_group_iterator = NodeGroups::iterator; using element_group_iterator = ElementGroups::iterator; using const_node_group_iterator = NodeGroups::const_iterator; using const_element_group_iterator = ElementGroups::const_iterator; #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(group_type, function, \ param_in, param_out) \ inline BOOST_PP_CAT(BOOST_PP_CAT(const_, group_type), _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) const { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ }; \ \ inline BOOST_PP_CAT(group_type, _iterator) \ BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) { \ return BOOST_PP_CAT(group_type, s).function(param_out); \ } #define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(group_type, function) \ AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION( \ group_type, function, BOOST_PP_EMPTY(), BOOST_PP_EMPTY()) AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, begin); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, end); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(element_group, find, const std::string & name, name); AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(node_group, find, const std::string & name, name); + public: decltype(auto) iterateNodeGroups() { return make_dereference_adaptor(make_values_adaptor(node_groups)); } decltype(auto) iterateNodeGroups() const { return make_dereference_adaptor(make_values_adaptor(node_groups)); } /* ------------------------------------------------------------------------ */ /* Clustering filter */ /* -------------------------------------------------------------------9+ ----- */ public: class ClusteringFilter { public: virtual bool operator()(const Element &) const { return true; } }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// create an empty node group NodeGroup & createNodeGroup(const std::string & group_name, bool replace_group = false); /// create a node group from another node group but filtered template NodeGroup & createFilteredNodeGroup(const std::string & group_name, const NodeGroup & node_group, T & filter); /// destroy a node group void destroyNodeGroup(const std::string & group_name); /// create an element group and the associated node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension = _all_dimensions, bool replace_group = false); /// create an element group from another element group but filtered template ElementGroup & createFilteredElementGroup(const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter); /// destroy an element group and the associated node group void destroyElementGroup(const std::string & group_name, bool destroy_node_group = false); /// destroy all element groups and the associated node groups void destroyAllElementGroups(bool destroy_node_groups = false); /// create a element group using an existing node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group); /// create groups based on values stored in a given mesh data template void createGroupsFromMeshData(const std::string & dataset_name); /// create boundaries group by a clustering algorithm \todo extend to parallel UInt createBoundaryGroupFromGeometry(); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, Mesh & mesh_facets, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); private: /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, const std::string & cluster_name_prefix, const ClusteringFilter & filter, Mesh & mesh_facets); public: /// Create an ElementGroup based on a NodeGroup void createElementGroupFromNodeGroup(const std::string & name, const std::string & node_group, UInt dimension = _all_dimensions); virtual void printself(std::ostream & stream, int indent = 0) const; /// this function insure that the group names are present on all processors /// /!\ it is a SMP call void synchronizeGroupNames(); -/// register an elemental field to the given group name (overloading for -/// ElementalPartionField) + /// register an elemental field to the given group name (overloading for + /// ElementalPartionField) template class dump_type> - dumper::Field * createElementalField( + std::shared_ptr createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// ElementalField) template class ret_type, template class, bool> class dump_type> - dumper::Field * createElementalField( + std::shared_ptr createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem = ElementTypeMap()); /// register an elemental field to the given group name (overloading for /// MaterialInternalField) template class dump_type> - dumper::Field * createElementalField(const ElementTypeMapArray & field, - const std::string & group_name, - UInt spatial_dimension, - const ElementKind & kind, - ElementTypeMap nb_data_per_elem); + std::shared_ptr + createElementalField(const ElementTypeMapArray & field, + const std::string & group_name, UInt spatial_dimension, + const ElementKind & kind, + ElementTypeMap nb_data_per_elem); template class ftype> - dumper::Field * createNodalField(const ftype * field, - const std::string & group_name, - UInt padding_size = 0); + std::shared_ptr + createNodalField(const ftype * field, + const std::string & group_name, UInt padding_size = 0); template class ftype> - dumper::Field * createStridedNodalField(const ftype * field, - const std::string & group_name, - UInt size, UInt stride, - UInt padding_size); + std::shared_ptr + createStridedNodalField(const ftype * field, + const std::string & group_name, UInt size, + UInt stride, UInt padding_size); protected: /// fill a buffer with all the group names void fillBufferWithGroupNames( CommunicationBufferTemplated & comm_buffer) const; /// take a buffer and create the missing groups localy void checkAndAddGroups(CommunicationBufferTemplated & buffer); /// register an elemental field to the given group name template - inline dumper::Field * + inline std::shared_ptr createElementalField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, const ElementTypeMap & nb_data_per_elem); /// register an elemental field to the given group name template - inline dumper::Field * + inline std::shared_ptr createElementalFilteredField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem); /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ElementGroups, element_groups, const ElementGroups &); const ElementGroup & getElementGroup(const std::string & name) const; const NodeGroup & getNodeGroup(const std::string & name) const; ElementGroup & getElementGroup(const std::string & name); NodeGroup & getNodeGroup(const std::string & name); UInt getNbElementGroups(UInt dimension = _all_dimensions) const; UInt getNbNodeGroups() { return node_groups.size(); }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id to create element and node groups ID id; /// memory_id to create element and node groups MemoryID memory_id; /// list of the node groups managed NodeGroups node_groups; /// list of the element groups managed ElementGroups element_groups; /// Mesh to which the element belongs const Mesh & mesh; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const GroupManager & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AKANTU_GROUP_MANAGER_HH__ */ diff --git a/src/mesh/group_manager_inline_impl.cc b/src/mesh/group_manager_inline_impl.cc index bfb9d0f6a..0b386bef1 100644 --- a/src/mesh/group_manager_inline_impl.cc +++ b/src/mesh/group_manager_inline_impl.cc @@ -1,205 +1,205 @@ /** * @file group_manager_inline_impl.cc * * @author Guillaume Anciaux * @author Dana Christen * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Sun Dec 03 2017 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * @section LICENSE * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dumper_field.hh" #include "element_group.hh" #include "element_type_map_filter.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_nodal_field.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template class dump_type> -dumper::Field * GroupManager::createElementalField( +std::shared_ptr GroupManager::createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const ElementTypeMapArray * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") return this->createElementalField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); else return this->createElementalFilteredField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); } /* -------------------------------------------------------------------------- */ template class T2, template class, bool> class dump_type> -dumper::Field * GroupManager::createElementalField( +std::shared_ptr GroupManager::createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const ElementTypeMapArray * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") return this->createElementalField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); else return this->createElementalFilteredField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); } /* -------------------------------------------------------------------------- */ template class dump_type> ///< type of InternalMaterialField -dumper::Field * -GroupManager::createElementalField(const ElementTypeMapArray & field, - const std::string & group_name, - UInt spatial_dimension, - const ElementKind & kind, - ElementTypeMap nb_data_per_elem) { +std::shared_ptr GroupManager::createElementalField( + const ElementTypeMapArray & field, const std::string & group_name, + UInt spatial_dimension, const ElementKind & kind, + ElementTypeMap nb_data_per_elem) { const ElementTypeMapArray * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") return this->createElementalField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); else return this->createElementalFilteredField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); } /* -------------------------------------------------------------------------- */ template -dumper::Field * GroupManager::createElementalField( +std::shared_ptr GroupManager::createElementalField( const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, const ElementTypeMap & nb_data_per_elem) { const field_type * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name != "all") throw; - dumper::Field * dumper = - new dump_type(field, spatial_dimension, _not_ghost, kind); + std::shared_ptr dumper = + std::make_shared(field, spatial_dimension, _not_ghost, kind); dumper->setNbDataPerElem(nb_data_per_elem); return dumper; } /* -------------------------------------------------------------------------- */ template -dumper::Field * GroupManager::createElementalFilteredField( +std::shared_ptr GroupManager::createElementalFilteredField( const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const field_type * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") throw; using T = typename field_type::type; ElementGroup & group = this->getElementGroup(group_name); UInt dim = group.getDimension(); if (dim != spatial_dimension) throw; const ElementTypeMapArray & elemental_filter = group.getElements(); auto * filtered = new ElementTypeMapArrayFilter(field, elemental_filter, nb_data_per_elem); - dumper::Field * dumper = new dump_type(*filtered, dim, _not_ghost, kind); + auto dumper = std::make_shared(*filtered, dim, _not_ghost, kind); dumper->setNbDataPerElem(nb_data_per_elem); return dumper; } /* -------------------------------------------------------------------------- */ template class ftype> -dumper::Field * GroupManager::createNodalField(const ftype * field, - const std::string & group_name, - UInt padding_size) { +std::shared_ptr +GroupManager::createNodalField(const ftype * field, + const std::string & group_name, + UInt padding_size) { if (field == nullptr) return nullptr; if (group_name == "all") { using DumpType = typename dumper::NodalField; - auto * dumper = new DumpType(*field, 0, 0, nullptr); + auto dumper = std::make_shared(*field, 0, 0, nullptr); dumper->setPadding(padding_size); return dumper; } else { ElementGroup & group = this->getElementGroup(group_name); const Array * nodal_filter = &(group.getNodes()); using DumpType = typename dumper::NodalField; - auto * dumper = new DumpType(*field, 0, 0, nodal_filter); + auto dumper = std::make_shared(*field, 0, 0, nodal_filter); dumper->setPadding(padding_size); return dumper; } return nullptr; } /* -------------------------------------------------------------------------- */ template class ftype> -dumper::Field * +std::shared_ptr GroupManager::createStridedNodalField(const ftype * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size) { if (field == NULL) return nullptr; if (group_name == "all") { using DumpType = typename dumper::NodalField; - auto * dumper = new DumpType(*field, size, stride, NULL); + auto dumper = std::make_shared(*field, size, stride, NULL); dumper->setPadding(padding_size); return dumper; } else { ElementGroup & group = this->getElementGroup(group_name); const Array * nodal_filter = &(group.getNodes()); using DumpType = typename dumper::NodalField; - auto * dumper = new DumpType(*field, size, stride, nodal_filter); + auto dumper = + std::make_shared(*field, size, stride, nodal_filter); dumper->setPadding(padding_size); return dumper; } return nullptr; } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc index c2d428059..71d99a6aa 100644 --- a/src/mesh/mesh.cc +++ b/src/mesh/mesh.cc @@ -1,575 +1,575 @@ /** * @file mesh.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief class handling 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 "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "group_manager_inline_impl.cc" #include "mesh.hh" #include "mesh_global_data_updater.hh" #include "mesh_io.hh" #include "mesh_iterators.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "element_synchronizer.hh" #include "facet_synchronizer.hh" #include "mesh_utils_distribution.hh" #include "node_synchronizer.hh" #include "periodic_node_synchronizer.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumper_field.hh" #include "dumper_internal_material_field.hh" #endif /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id, Communicator & communicator) : Memory(id, memory_id), GroupManager(*this, id + ":group_manager", memory_id), MeshData("mesh_data", id, memory_id), connectivities("connectivities", id, memory_id), ghosts_counters("ghosts_counters", id, memory_id), normals("normals", id, memory_id), spatial_dimension(spatial_dimension), size(spatial_dimension, 0.), bbox(spatial_dimension), bbox_local(spatial_dimension), communicator(&communicator) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, communicator) { AKANTU_DEBUG_IN(); this->nodes = std::make_shared>(0, spatial_dimension, id + ":coordinates"); this->nodes_flags = std::make_shared>(0, 1, NodeFlag::_normal, id + ":nodes_flags"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id, memory_id) {} /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id, const MemoryID & memory_id) : Mesh(spatial_dimension, id, memory_id, Communicator::getStaticCommunicator()) { this->nodes = nodes; this->nb_global_nodes = this->nodes->size(); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getBarycenters(Array & barycenter, const ElementType & type, const GhostType & ghost_type) const { barycenter.resize(getNbElement(type, ghost_type)); for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { getBarycenter(Element{type, UInt(std::get<0>(data)), ghost_type}, std::get<1>(data)); } } /* -------------------------------------------------------------------------- */ Mesh & Mesh::initMeshFacets(const ID & id) { AKANTU_DEBUG_IN(); if (mesh_facets) { AKANTU_DEBUG_OUT(); return *mesh_facets; } mesh_facets = std::make_unique(spatial_dimension, this->nodes, getID() + ":" + id, getMemoryID()); mesh_facets->mesh_parent = this; mesh_facets->is_mesh_facets = true; mesh_facets->nodes_flags = this->nodes_flags; mesh_facets->nodes_global_ids = this->nodes_global_ids; MeshUtils::buildAllFacets(*this, *mesh_facets, 0); if (mesh.isDistributed()) { mesh_facets->is_distributed = true; mesh_facets->element_synchronizer = std::make_unique( *mesh_facets, mesh.getElementSynchronizer()); } /// transfers the the mesh physical names to the mesh facets if (not this->hasData("physical_names")) { AKANTU_DEBUG_OUT(); return *mesh_facets; } if (not mesh_facets->hasData("physical_names")) { mesh_facets->registerElementalData("physical_names"); } auto & mesh_phys_data = this->getData("physical_names"); auto & phys_data = mesh_facets->getData("physical_names"); phys_data.initialize(*mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); ElementTypeMapArray barycenters(getID(), "temporary_barycenters"); barycenters.initialize(*mesh_facets, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); for (auto && ghost_type : ghost_types) { for (auto && type : barycenters.elementTypes(spatial_dimension - 1, ghost_type)) { mesh_facets->getBarycenters(barycenters(type, ghost_type), type, ghost_type); } } for_each_element( mesh, [&](auto && element) { Vector barycenter(spatial_dimension); mesh.getBarycenter(element, barycenter); auto norm_barycenter = barycenter.norm(); auto tolerance = Math::getTolerance(); if (norm_barycenter > tolerance) tolerance *= norm_barycenter; const auto & element_to_facet = mesh_facets->getElementToSubelement( element.type, element.ghost_type); Vector barycenter_facet(spatial_dimension); auto range = enumerate(make_view( barycenters(element.type, element.ghost_type), spatial_dimension)); #ifndef AKANTU_NDEBUG auto min_dist = std::numeric_limits::max(); #endif // this is a spacial search coded the most inefficient way. auto facet = std::find_if(range.begin(), range.end(), [&](auto && data) { auto facet = std::get<0>(data); if (element_to_facet(facet)[1] == ElementNull) return false; auto norm_distance = barycenter.distance(std::get<1>(data)); #ifndef AKANTU_NDEBUG min_dist = std::min(min_dist, norm_distance); #endif return (norm_distance < tolerance); }); if (facet == range.end()) { AKANTU_DEBUG_INFO("The element " << element << " did not find its associated facet in the " "mesh_facets! Try to decrease math tolerance. " "The closest element was at a distance of " << min_dist); return; } // set physical name phys_data(Element{element.type, UInt(std::get<0>(*facet)), element.ghost_type}) = mesh_phys_data(element); }, _spatial_dimension = spatial_dimension - 1); mesh_facets->createGroupsFromMeshData("physical_names"); AKANTU_DEBUG_OUT(); return *mesh_facets; } /* -------------------------------------------------------------------------- */ void Mesh::defineMeshParent(const Mesh & mesh) { AKANTU_DEBUG_IN(); this->mesh_parent = &mesh; this->is_mesh_facets = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() = default; /* -------------------------------------------------------------------------- */ void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) { AKANTU_DEBUG_ASSERT(not is_distributed, "You cannot read a mesh that is already distributed"); MeshIO mesh_io; mesh_io.read(filename, *this, mesh_io_type); auto types = this->elementTypes(spatial_dimension, _not_ghost, _ek_not_defined); auto it = types.begin(); auto last = types.end(); if (it == last) AKANTU_DEBUG_WARNING( "The mesh contained in the file " << filename << " does not seem to be of the good dimension." << " No element of dimension " << spatial_dimension << " where read."); this->makeReady(); } /* -------------------------------------------------------------------------- */ void Mesh::write(const std::string & filename, const MeshIOType & mesh_io_type) { MeshIO mesh_io; mesh_io.write(filename, *this, mesh_io_type); } /* -------------------------------------------------------------------------- */ void Mesh::makeReady() { this->nb_global_nodes = this->nodes->size(); this->computeBoundingBox(); this->nodes_flags->resize(nodes->size(), NodeFlag::_normal); this->nodes_to_elements.resize(nodes->size()); for (auto & node_set : nodes_to_elements) { node_set = std::make_unique>(); } } /* -------------------------------------------------------------------------- */ void Mesh::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Mesh [" << std::endl; stream << space << " + id : " << getID() << std::endl; stream << space << " + spatial dimension : " << this->spatial_dimension << std::endl; stream << space << " + nodes [" << std::endl; nodes->printself(stream, indent + 2); stream << space << " + connectivities [" << std::endl; connectivities.printself(stream, indent + 2); stream << space << " ]" << std::endl; GroupManager::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Mesh::computeBoundingBox() { AKANTU_DEBUG_IN(); bbox_local.reset(); for (auto & pos : make_view(*nodes, spatial_dimension)) { // if(!isPureGhostNode(i)) bbox_local += pos; } if (this->is_distributed) { bbox = bbox_local.allSum(*communicator); } else { bbox = bbox_local; } size = bbox.size(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initNormals() { normals.initialize(*this, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _element_kind = _ek_not_defined); } /* -------------------------------------------------------------------------- */ void Mesh::getGlobalConnectivity( ElementTypeMapArray & global_connectivity) { AKANTU_DEBUG_IN(); for (auto && ghost_type : ghost_types) { for (auto type : global_connectivity.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_not_defined, _ghost_type = ghost_type)) { if (not connectivities.exists(type, ghost_type)) continue; auto & local_conn = connectivities(type, ghost_type); auto & g_connectivity = global_connectivity(type, ghost_type); UInt nb_nodes = local_conn.size() * local_conn.getNbComponent(); std::transform(local_conn.begin_reinterpret(nb_nodes), local_conn.end_reinterpret(nb_nodes), g_connectivity.begin_reinterpret(nb_nodes), [&](UInt l) -> UInt { return this->getNodeGlobalId(l); }); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Mesh::getGroupDumper(const std::string & dumper_name, const std::string & group_name) { if (group_name == "all") return this->getDumper(dumper_name); else return element_groups[group_name]->getDumper(dumper_name); } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & arrays, const ElementKind & element_kind) { ElementTypeMap nb_data_per_elem; for (auto type : elementTypes(spatial_dimension, _not_ghost, element_kind)) { UInt nb_elements = this->getNbElement(type); auto & array = arrays(type); nb_data_per_elem(type) = array.getNbComponent() * array.size(); nb_data_per_elem(type) /= nb_elements; } return nb_data_per_elem; } /* -------------------------------------------------------------------------- */ template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); template ElementTypeMap Mesh::getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER template -dumper::Field * +std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind) { - dumper::Field * field = nullptr; + std::shared_ptr field; ElementTypeMapArray * internal = nullptr; try { internal = &(this->getData(field_id)); } catch (...) { return nullptr; } ElementTypeMap nb_data_per_elem = this->getNbDataPerElem(*internal, element_kind); field = this->createElementalField( *internal, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); return field; } -template dumper::Field * +template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); -template dumper::Field * +template std::shared_ptr Mesh::createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, const ElementKind & element_kind); #endif /* -------------------------------------------------------------------------- */ void Mesh::distributeImpl( Communicator & communicator, std::function edge_weight_function, std::function vertex_weight_function) { AKANTU_DEBUG_ASSERT(is_distributed == false, "This mesh is already distribute"); this->communicator = &communicator; this->element_synchronizer = std::make_unique( *this, this->getID() + ":element_synchronizer", this->getMemoryID(), true); this->node_synchronizer = std::make_unique( *this, this->getID() + ":node_synchronizer", this->getMemoryID(), true); Int psize = this->communicator->getNbProc(); if (psize > 1) { #ifdef AKANTU_USE_SCOTCH Int prank = this->communicator->whoAmI(); if (prank == 0) { MeshPartitionScotch partition(*this, spatial_dimension); partition.partitionate(psize, edge_weight_function, vertex_weight_function); MeshUtilsDistribution::distributeMeshCentralized(*this, 0, partition); } else { MeshUtilsDistribution::distributeMeshCentralized(*this, 0); } #else if (psize > 1) { AKANTU_ERROR("Cannot distribute a mesh without a partitioning tool"); } #endif } // if (psize > 1) this->is_distributed = true; this->computeBoundingBox(); } /* -------------------------------------------------------------------------- */ void Mesh::getAssociatedElements(const Array & node_list, Array & elements) { for (const auto & node : node_list) for (const auto & element : *nodes_to_elements[node]) elements.push_back(element); } /* -------------------------------------------------------------------------- */ void Mesh::fillNodesToElements() { Element e; UInt nb_nodes = nodes->size(); for (UInt n = 0; n < nb_nodes; ++n) { if (this->nodes_to_elements[n]) this->nodes_to_elements[n]->clear(); else this->nodes_to_elements[n] = std::make_unique>(); } for (auto ghost_type : ghost_types) { e.ghost_type = ghost_type; for (const auto & type : elementTypes(spatial_dimension, ghost_type, _ek_not_defined)) { e.type = type; UInt nb_element = this->getNbElement(type, ghost_type); Array::const_iterator> conn_it = connectivities(type, ghost_type) .begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; const Vector & conn = *conn_it; for (UInt n = 0; n < conn.size(); ++n) nodes_to_elements[conn(n)]->insert(e); } } } } /* -------------------------------------------------------------------------- */ std::tuple Mesh::updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) { if (global_data_updater) return this->global_data_updater->updateData(nodes_event, elements_event); else { return std::make_tuple(nodes_event.getList().size(), elements_event.getList().size()); } } /* -------------------------------------------------------------------------- */ void Mesh::registerGlobalDataUpdater( std::unique_ptr && global_data_updater) { this->global_data_updater = std::move(global_data_updater); } /* -------------------------------------------------------------------------- */ void Mesh::eraseElements(const Array & elements) { ElementTypeMap last_element; RemovedElementsEvent event(*this); auto & remove_list = event.getList(); auto & new_numbering = event.getNewNumbering(); for (auto && el : elements) { if (el.ghost_type != _not_ghost) { auto & count = ghosts_counters(el); --count; if (count > 0) continue; } remove_list.push_back(el); if (not last_element.exists(el.type, el.ghost_type)) { UInt nb_element = mesh.getNbElement(el.type, el.ghost_type); last_element(nb_element - 1, el.type, el.ghost_type); auto & numbering = new_numbering.alloc(nb_element, 1, el.type, el.ghost_type); for (auto && pair : enumerate(numbering)) { std::get<1>(pair) = std::get<0>(pair); } } UInt & pos = last_element(el.type, el.ghost_type); auto & numbering = new_numbering(el.type, el.ghost_type); numbering(el.element) = UInt(-1); numbering(pos) = el.element; --pos; } this->sendEvent(event); } } // namespace akantu diff --git a/src/mesh/mesh.hh b/src/mesh/mesh.hh index 8da337278..ce26878b3 100644 --- a/src/mesh/mesh.hh +++ b/src/mesh/mesh.hh @@ -1,695 +1,696 @@ /** * @file mesh.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Feb 19 2018 * * @brief the class representing the 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_bbox.hh" #include "aka_event_handler_manager.hh" #include "aka_memory.hh" #include "communicator.hh" #include "dumpable.hh" #include "element.hh" #include "element_class.hh" #include "element_type_map.hh" #include "group_manager.hh" #include "mesh_data.hh" #include "mesh_events.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { class ElementSynchronizer; class NodeSynchronizer; class PeriodicNodeSynchronizer; class MeshGlobalDataUpdater; } // namespace akantu namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(communicator); DECLARE_NAMED_ARGUMENT(edge_weight_function); DECLARE_NAMED_ARGUMENT(vertex_weight_function); } // namespace /* -------------------------------------------------------------------------- */ /* Mesh */ /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes * Array, and the connectivity. The connectivity are stored in by element * types. * * In order to loop on all element you have to loop on all types like this : * @code{.cpp} for(auto & type : mesh.elementTypes()) { UInt nb_element = mesh.getNbElement(type); const Array & conn = mesh.getConnectivity(type); for(UInt e = 0; e < nb_element; ++e) { ... } } or for_each_element(mesh, [](Element & element) { std::cout << element << std::endl }); @endcode */ class Mesh : protected Memory, public EventHandlerManager, public GroupManager, public MeshData, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: /// default constructor used for chaining, the last parameter is just to /// differentiate constructors Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id, Communicator & communicator); public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const ID & id = "mesh", const MemoryID & memory_id = 0); /// mesh not distributed and not using the default communicator Mesh(UInt spatial_dimension, Communicator & communicator, const ID & id = "mesh", const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, const std::shared_ptr> & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); ~Mesh() override; /// read the mesh from a file void read(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); /// write the mesh to a file void write(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); protected: void makeReady(); private: /// initialize the connectivity to NULL and other stuff void init(); /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /* ------------------------------------------------------------------------ */ /* Distributed memory methods and accessors */ /* ------------------------------------------------------------------------ */ public: protected: /// patitionate the mesh among the processors involved in their computation virtual void distributeImpl( Communicator & communicator, std::function edge_weight_function, std::function vertex_weight_function); public: /// with the arguments to pass to the partitionner template std::enable_if_t::value> distribute(pack &&... _pack) { distributeImpl( OPTIONAL_NAMED_ARG(communicator, Communicator::getStaticCommunicator()), OPTIONAL_NAMED_ARG(edge_weight_function, [](auto &&, auto &&) { return 1; }), OPTIONAL_NAMED_ARG(vertex_weight_function, [](auto &&) { return 1; })); } /// defines is the mesh is distributed or not inline bool isDistributed() const { return this->is_distributed; } /* ------------------------------------------------------------------------ */ /* Periodicity methods and accessors */ /* ------------------------------------------------------------------------ */ public: /// set the periodicity in a given direction void makePeriodic(const SpatialDirection & direction); void makePeriodic(const SpatialDirection & direction, const ID & list_1, const ID & list_2); protected: void makePeriodic(const SpatialDirection & direction, const Array & list_1, const Array & list_2); /// Removes the face that the mesh is periodic void wipePeriodicInfo(); inline void addPeriodicSlave(UInt slave, UInt master); template void synchronizePeriodicSlaveDataWithMaster(Array & data); // update the periodic synchronizer (creates it if it does not exists) void updatePeriodicSynchronizer(); public: /// defines if the mesh is periodic or not inline bool isPeriodic() const { return (this->is_periodic != 0); } inline bool isPeriodic(const SpatialDirection & direction) const { return ((this->is_periodic & (1 << direction)) != 0); } class PeriodicSlaves; /// get the master node for a given slave nodes, except if node not a slave inline UInt getPeriodicMaster(UInt slave) const; /// get an iterable list of slaves for a given master node inline decltype(auto) getPeriodicSlaves(UInt master) const; /* ------------------------------------------------------------------------ */ /* General Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Array & nodal_values, T * elemental_values, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const; // /// extract coordinates of nodes from a reversed element // inline void extractNodalCoordinatesFromPBCElement(Real * local_coords, // UInt * connectivity, // UInt n_nodes); /// add a Array of connectivity for the type . inline void addConnectivityType(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ template inline void sendEvent(Event & event) { // if(event.getList().size() != 0) EventHandlerManager::sendEvent(event); } /// prepare the event to remove the elements listed void eraseElements(const Array & elements); /* ------------------------------------------------------------------------ */ template inline void removeNodesFromArray(Array & vect, const Array & new_numbering); /// initialize normals void initNormals(); /// init facets' mesh Mesh & initMeshFacets(const ID & id = "mesh_facets"); /// define parent mesh void defineMeshParent(const Mesh & mesh); /// get global connectivity array void getGlobalConnectivity(ElementTypeMapArray & global_connectivity); public: void getAssociatedElements(const Array & node_list, Array & elements); private: /// fills the nodes_to_elements structure void fillNodesToElements(); /// update the global ids, nodes type, ... std::tuple updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event); void registerGlobalDataUpdater( std::unique_ptr && global_data_updater); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the id of the mesh AKANTU_GET_MACRO(ID, Memory::id, const ID &); /// get the id of the mesh AKANTU_GET_MACRO(MemoryID, Memory::memory_id, const MemoryID &); /// get the spatial dimension of the mesh = number of component of the /// coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Array aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Array &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Array &); /// get the normals for the elements AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, Real); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->size(), UInt); /// get the Array of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Array &); // AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array // &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global id of a node inline UInt getNodeLocalId(UInt global_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Array AKANTU_GET_MACRO(NodesFlags, *nodes_flags, const Array &); protected: AKANTU_GET_MACRO_NOT_CONST(NodesFlags, *nodes_flags, Array &); public: inline NodeFlag getNodeFlag(UInt local_id) const; inline Int getNodePrank(UInt local_id) const; /// say if a node is a pure ghost node inline bool isPureGhostNode(UInt n) const; /// say if a node is pur local or master node inline bool isLocalOrMasterNode(UInt n) const; inline bool isLocalNode(UInt n) const; inline bool isMasterNode(UInt n) const; inline bool isSlaveNode(UInt n) const; inline bool isPeriodicSlave(UInt n) const; inline bool isPeriodicMaster(UInt n) const; const Vector & getLowerBounds() const { return bbox.getLowerBounds(); } const Vector & getUpperBounds() const { return bbox.getUpperBounds(); } AKANTU_GET_MACRO(BBox, bbox, const BBox &); const Vector & getLocalLowerBounds() const { return bbox_local.getLowerBounds(); } const Vector & getLocalUpperBounds() const { return bbox_local.getUpperBounds(); } AKANTU_GET_MACRO(LocalBBox, bbox_local, const BBox &); /// get the connectivity Array for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); AKANTU_GET_MACRO(Connectivities, connectivities, const ElementTypeMapArray &); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// get the number of element for a given ghost_type and a given dimension inline UInt getNbElement(const UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_not_defined) const; /// compute the barycenter of a given element inline void getBarycenter(const Element & element, Vector & barycenter) const; void getBarycenters(Array & barycenter, const ElementType & type, const GhostType & ghost_type) const; /// get the element connected to a subelement (element of lower dimension) const auto & getElementToSubelement() const; /// get the element connected to a subelement const auto & getElementToSubelement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get the element connected to a subelement auto & getElementToSubelement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the elements connected to a subelement const auto & getElementToSubelement(const Element & element) const; /// get the subelement (element of lower dimension) connected to a element const auto & getSubelementToElement() const; /// get the subelement connected to an element const auto & getSubelementToElement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get the subelement connected to an element auto & getSubelementToElement(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the subelement (element of lower dimension) connected to a element VectorProxy getSubelementToElement(const Element & element) const; /// get connectivity of a given element inline VectorProxy getConnectivity(const Element & element) const; inline Vector getConnectivityWithPeriodicity(const Element & element) const; protected: inline auto & getElementToSubelement(const Element & element); inline VectorProxy getSubelementToElement(const Element & element); inline VectorProxy getConnectivity(const Element & element); public: /// get a name field associated to the mesh template inline const Array & getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get a name field associated to the mesh template inline Array & getData(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get a name field associated to the mesh template inline const ElementTypeMapArray & getData(const ID & data_name) const; /// get a name field associated to the mesh template inline ElementTypeMapArray & getData(const ID & data_name); template ElementTypeMap getNbDataPerElem(ElementTypeMapArray & array, const ElementKind & element_kind); template - dumper::Field * createFieldFromAttachedData(const std::string & field_id, - const std::string & group_name, - const ElementKind & element_kind); + std::shared_ptr + createFieldFromAttachedData(const std::string & field_id, + const std::string & group_name, + const ElementKind & element_kind); /// templated getter returning the pointer to data in MeshData (modifiable) template inline Array & getDataPointer(const std::string & data_name, const ElementType & el_type, const GhostType & ghost_type = _not_ghost, UInt nb_component = 1, bool size_to_nb_element = true, bool resize_with_parent = false); template inline Array & getDataPointer(const ID & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent, const T & defaul_); /// Facets mesh accessor inline const Mesh & getMeshFacets() const; inline Mesh & getMeshFacets(); /// Parent mesh accessor inline const Mesh & getMeshParent() const; inline bool isMeshFacets() const { return this->is_mesh_facets; } /// return the dumper from a group and and a dumper name DumperIOHelper & getGroupDumper(const std::string & dumper_name, const std::string & group_name); /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get the kind of the element type static inline ElementKind getKind(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type, UInt t); /// get local connectivity of a facet for a given facet type static inline auto getFacetLocalConnectivity(const ElementType & type, UInt t = 0); /// get connectivity of facets for a given element inline auto getFacetConnectivity(const Element & element, UInt t = 0) const; /// get the number of type of the surface element associated to a given /// element type static inline UInt getNbFacetTypes(const ElementType & type, UInt t = 0); /// get the type of the surface element associated to a given element static inline constexpr auto getFacetType(const ElementType & type, UInt t = 0); /// get all the type of the surface element associated to a given element static inline constexpr auto getAllFacetTypes(const ElementType & type); /// get the number of nodes in the given element list static inline UInt getNbNodesPerElementList(const Array & elements); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ using type_iterator [[deprecated]] = ElementTypeMapArray::type_iterator; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; template ElementTypesIteratorHelper elementTypes(pack &&... _pack) const; [[deprecated("Use elementTypes instead")]] inline decltype(auto) firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.elementTypes(dim, ghost_type, kind).begin(); } [[deprecated("Use elementTypes instead")]] inline decltype(auto) lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.elementTypes(dim, ghost_type, kind).end(); } AKANTU_GET_MACRO(ElementSynchronizer, *element_synchronizer, const ElementSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(ElementSynchronizer, *element_synchronizer, ElementSynchronizer &); AKANTU_GET_MACRO(NodeSynchronizer, *node_synchronizer, const NodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(NodeSynchronizer, *node_synchronizer, NodeSynchronizer &); AKANTU_GET_MACRO(PeriodicNodeSynchronizer, *periodic_node_synchronizer, const PeriodicNodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(PeriodicNodeSynchronizer, *periodic_node_synchronizer, PeriodicNodeSynchronizer &); -// AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator -// &); + // AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator + // &); AKANTU_GET_MACRO(Communicator, *communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, auto &); AKANTU_GET_MACRO(PeriodicMasterSlaves, periodic_master_slave, const auto &); /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshAccessor; friend class MeshUtils; AKANTU_GET_MACRO(NodesPointer, *nodes, Array &); /// get a pointer to the nodes_global_ids Array and create it if /// necessary inline Array & getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Array and create it if necessary inline Array & getNodesFlagsPointer(); /// get a pointer to the connectivity Array for the given type and create it /// if necessary inline Array & getConnectivityPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get the ghost element counter inline Array & getGhostsCounters(const ElementType & type, const GhostType & ghost_type = _ghost) { AKANTU_DEBUG_ASSERT(ghost_type != _not_ghost, "No ghost counter for _not_ghost elements"); return ghosts_counters(type, ghost_type); } /// get a pointer to the element_to_subelement Array for the given type and /// create it if necessary inline Array> & getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Array for the given type and /// create it if necessary inline Array & getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// array of the nodes coordinates std::shared_ptr> nodes; /// global node ids std::shared_ptr> nodes_global_ids; /// node flags (shared/periodic/...) std::shared_ptr> nodes_flags; /// processor handling the node when not local or master std::unordered_map nodes_prank; /// global number of nodes; UInt nb_global_nodes{0}; /// all class of elements present in this mesh (for heterogenous meshes) ElementTypeMapArray connectivities; /// count the references on ghost elements ElementTypeMapArray ghosts_counters; /// map to normals for all class of elements present in this mesh ElementTypeMapArray normals; /// the spatial dimension of this mesh UInt spatial_dimension{0}; /// size covered by the mesh on each direction Vector size; /// global bounding box BBox bbox; /// local bounding box BBox bbox_local; /// Extra data loaded from the mesh file // MeshData mesh_data; /// facets' mesh std::unique_ptr mesh_facets; /// parent mesh (this is set for mesh_facets meshes) const Mesh * mesh_parent{nullptr}; /// defines if current mesh is mesh_facets or not bool is_mesh_facets{false}; /// defines if the mesh is centralized or distributed bool is_distributed{false}; /// defines if the mesh is periodic bool is_periodic{false}; /// Communicator on which mesh is distributed Communicator * communicator; /// Element synchronizer std::unique_ptr element_synchronizer; /// Node synchronizer std::unique_ptr node_synchronizer; /// Node synchronizer for periodic nodes std::unique_ptr periodic_node_synchronizer; using NodesToElements = std::vector>>; /// class to update global data using external knowledge std::unique_ptr global_data_updater; /// This info is stored to simplify the dynamic changes NodesToElements nodes_to_elements; /// periodicity local info std::unordered_map periodic_slave_master; std::unordered_multimap periodic_master_slave; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ #include "element_type_map_tmpl.hh" #include "mesh_inline_impl.cc" #endif /* __AKANTU_MESH_HH__ */ diff --git a/src/mesh/node_group.cc b/src/mesh/node_group.cc index 5ece1d1d0..6b3f0a02f 100644 --- a/src/mesh/node_group.cc +++ b/src/mesh/node_group.cc @@ -1,98 +1,98 @@ /** * @file node_group.cc * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Feb 01 2018 * * @brief Implementation of the node group * * @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 "node_group.hh" #include "dumpable.hh" #include "dumpable_inline_impl.hh" #include "mesh.hh" #if defined(AKANTU_USE_IOHELPER) #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NodeGroup::NodeGroup(const std::string & name, const Mesh & mesh, const std::string & id, const MemoryID & memory_id) : Memory(id, memory_id), name(name), node_group(alloc(std::string(this->id + ":nodes"), 0, 1)) { #if defined(AKANTU_USE_IOHELPER) this->registerDumper("paraview_" + name, name, true); - this->getDumper().registerField( - "positions", new dumper::NodalField(mesh.getNodes(), 0, 0, - &this->getNodes())); + auto field = std::make_shared>( + mesh.getNodes(), 0, 0, &this->getNodes()); + this->getDumper().registerField("positions", field); #endif } /* -------------------------------------------------------------------------- */ NodeGroup::~NodeGroup() = default; /* -------------------------------------------------------------------------- */ void NodeGroup::empty() { node_group.resize(0); } /* -------------------------------------------------------------------------- */ void NodeGroup::optimize() { std::sort(node_group.begin(), node_group.end()); Array::iterator<> end = std::unique(node_group.begin(), node_group.end()); node_group.resize(end - node_group.begin()); } /* -------------------------------------------------------------------------- */ void NodeGroup::append(const NodeGroup & other_group) { AKANTU_DEBUG_IN(); UInt nb_nodes = node_group.size(); /// append new nodes to current list node_group.resize(nb_nodes + other_group.node_group.size()); std::copy(other_group.node_group.begin(), other_group.node_group.end(), node_group.begin() + nb_nodes); optimize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NodeGroup::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NodeGroup [" << std::endl; stream << space << " + name: " << name << std::endl; node_group.printself(stream, indent + 1); stream << space << "]" << std::endl; } -} // akantu +} // namespace akantu diff --git a/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index 1a1838901..7801f9db9 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,953 +1,961 @@ /** * @file heat_transfer_model.cc * * @author Guillaume Anciaux * @author Lucas Frerot * @author Emil Gallyamov * @author David Simon Kammer * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Tue Feb 20 2018 * * @brief Implementation of HeatTransferModel class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_model.hh" #include "dumpable_inline_impl.hh" #include "element_synchronizer.hh" #include "fe_engine_template.hh" #include "generalized_trapezoidal.hh" #include "group_manager_inline_impl.cc" #include "integrator_gauss.hh" #include "mesh.hh" #include "parser.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { namespace heat_transfer { namespace details { class ComputeRhoFunctor { public: ComputeRhoFunctor(const HeatTransferModel & model) : model(model){}; void operator()(Matrix & rho, const Element &) const { - rho.set(model.getCapacity()*model.getDensity()); + rho.set(model.getCapacity() * model.getDensity()); } private: const HeatTransferModel & model; }; - } -} + } // namespace details +} // namespace heat_transfer /* -------------------------------------------------------------------------- */ HeatTransferModel::HeatTransferModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, ModelType::_heat_transfer_model, dim, id, memory_id), temperature_gradient("temperature_gradient", id), temperature_on_qpoints("temperature_on_qpoints", id), conductivity_on_qpoints("conductivity_on_qpoints", id), k_gradt_on_qpoints("k_gradt_on_qpoints", id) { AKANTU_DEBUG_IN(); conductivity = Matrix(this->spatial_dimension, this->spatial_dimension); this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, _gst_htm_temperature); this->registerSynchronizer(synchronizer, _gst_htm_gradient_temperature); } registerFEEngineObject(id + ":fem", mesh, spatial_dimension); #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("heat_transfer", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif this->registerParam("conductivity", conductivity, _pat_parsmod); this->registerParam("conductivity_variation", conductivity_variation, 0., _pat_parsmod); this->registerParam("temperature_reference", T_ref, 0., _pat_parsmod); this->registerParam("capacity", capacity, _pat_parsmod); this->registerParam("density", density, _pat_parsmod); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initModel() { auto & fem = this->getFEEngine(); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); temperature_on_qpoints.initialize(fem, _nb_component = 1); temperature_gradient.initialize(fem, _nb_component = spatial_dimension); - conductivity_on_qpoints.initialize( - fem, _nb_component = spatial_dimension * spatial_dimension); + conductivity_on_qpoints.initialize(fem, _nb_component = spatial_dimension * + spatial_dimension); k_gradt_on_qpoints.initialize(fem, _nb_component = spatial_dimension); } /* -------------------------------------------------------------------------- */ FEEngine & HeatTransferModel::getFEEngineBoundary(const ID & name) { return dynamic_cast(getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::allocNodalField(Array *& array, const ID & name) { if (array == nullptr) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":" << name; array = &(alloc(sstr_disp.str(), nb_nodes, 1, T())); } } /* -------------------------------------------------------------------------- */ HeatTransferModel::~HeatTransferModel() = default; /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor compute_rho(*this); - for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { + for (auto & type : + mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldLumped(compute_rho, "M", "temperature", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType HeatTransferModel::getMatrixType(const ID & matrix_id) { if (matrix_id == "K" or matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleConductivityMatrix(); } else if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacity(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacityLumped(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleResidual() { AKANTU_DEBUG_IN(); this->assembleInternalHeatRate(); this->getDOFManager().assembleToResidual("temperature", *this->external_heat_rate, 1); this->getDOFManager().assembleToResidual("temperature", *this->internal_heat_rate, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::predictor() { ++temperature_release; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped() { AKANTU_DEBUG_IN(); if (!this->getDOFManager().hasLumpedMatrix("M")) { this->getDOFManager().getNewLumpedMatrix("M"); } this->getDOFManager().clearLumpedMatrix("M"); assembleCapacityLumped(_not_ghost); assembleCapacityLumped(_ghost); need_to_reassemble_capacity_lumped = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->temperature, "temperature"); this->allocNodalField(this->external_heat_rate, "external_heat_rate"); this->allocNodalField(this->internal_heat_rate, "internal_heat_rate"); this->allocNodalField(this->blocked_dofs, "blocked_dofs"); if (!dof_manager.hasDOFs("temperature")) { dof_manager.registerDOFs("temperature", *this->temperature, _dst_nodal); dof_manager.registerBlockedDOFs("temperature", *this->blocked_dofs); } if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(this->temperature_rate, "temperature_rate"); if (!dof_manager.hasDOFsDerivatives("temperature", 1)) { dof_manager.registerDOFsDerivative("temperature", 1, *this->temperature_rate); } } } /* -------------------------------------------------------------------------- */ std::tuple HeatTransferModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", _tsst_dynamic_lumped); } case _static: { return std::make_tuple("static", _tsst_static); } case _implicit_dynamic: { return std::make_tuple("implicit", _tsst_dynamic); } default: return std::make_tuple("unknown", _tsst_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions HeatTransferModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic_lumped: { options.non_linear_solver_type = _nls_lumped; options.integration_scheme_type["temperature"] = _ist_forward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["temperature"] = _ist_pseudo_time; options.solution_type["temperature"] = IntegrationScheme::_not_defined; break; } case _tsst_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["temperature"] = _ist_forward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; } else { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["temperature"] = _ist_backward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleConductivityMatrix() { AKANTU_DEBUG_IN(); this->computeConductivityOnQuadPoints(_not_ghost); if (conductivity_release[_not_ghost] == conductivity_matrix_release) return; if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().clearMatrix("K"); switch (mesh.getSpatialDimension()) { case 1: this->assembleConductivityMatrix<1>(_not_ghost); break; case 2: this->assembleConductivityMatrix<2>(_not_ghost); break; case 3: this->assembleConductivityMatrix<3>(_not_ghost); break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::assembleConductivityMatrix( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine(); - for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { + for (auto && type : + mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto nb_element = mesh.getNbElement(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); fem.computeBtDB(conductivity_on_qpoints(type, ghost_type), *bt_d_b, 2, type, ghost_type); /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto K_e = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_e"); fem.integrate(*bt_d_b, *K_e, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "temperature", *K_e, type, ghost_type, _symmetric); } conductivity_matrix_release = conductivity_release[ghost_type]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeConductivityOnQuadPoints( const GhostType & ghost_type) { // if already computed once check if need to compute if (not initial_conductivity[ghost_type]) { // if temperature did not change, condictivity will not vary if (temperature_release == conductivity_release[ghost_type]) return; // if conductivity_variation is 0 no need to recompute if (conductivity_variation == 0.) return; } - for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { + for (auto & type : + mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & temperature_interpolated = temperature_on_qpoints(type, ghost_type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type, ghost_type); auto & cond = conductivity_on_qpoints(type, ghost_type); for (auto && tuple : zip(make_view(cond, spatial_dimension, spatial_dimension), temperature_interpolated)) { auto & C = std::get<0>(tuple); auto & T = std::get<1>(tuple); C = conductivity; Matrix variation(spatial_dimension, spatial_dimension, conductivity_variation * (T - T_ref)); // @TODO: Guillaume are you sure ? why due you compute variation then ? C += conductivity_variation; } } conductivity_release[ghost_type] = temperature_release; initial_conductivity[ghost_type] = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeKgradT(const GhostType & ghost_type) { computeConductivityOnQuadPoints(ghost_type); - for (auto & type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { + for (auto & type : + mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & gradient = temperature_gradient(type, ghost_type); this->getFEEngine().gradientOnIntegrationPoints(*temperature, gradient, 1, type, ghost_type); for (auto && values : zip(make_view(conductivity_on_qpoints(type, ghost_type), spatial_dimension, spatial_dimension), make_view(gradient, spatial_dimension), make_view(k_gradt_on_qpoints(type, ghost_type), spatial_dimension))) { const auto & C = std::get<0>(values); const auto & BT = std::get<1>(values); auto & k_BT = std::get<2>(values); k_BT.mul(C, BT); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleInternalHeatRate() { AKANTU_DEBUG_IN(); this->internal_heat_rate->clear(); this->synchronize(_gst_htm_temperature); auto & fem = this->getFEEngine(); for (auto ghost_type : ghost_types) { // compute k \grad T computeKgradT(ghost_type); - for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { + for (auto type : + mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto & k_gradt_on_qpoints_vect = k_gradt_on_qpoints(type, ghost_type); UInt nb_quad_points = k_gradt_on_qpoints_vect.size(); Array bt_k_gT(nb_quad_points, nb_nodes_per_element); fem.computeBtD(k_gradt_on_qpoints_vect, bt_k_gT, type, ghost_type); UInt nb_elements = mesh.getNbElement(type, ghost_type); Array int_bt_k_gT(nb_elements, nb_nodes_per_element); fem.integrate(bt_k_gT, int_bt_k_gT, nb_nodes_per_element, type, ghost_type); this->getDOFManager().assembleElementalArrayLocalArray( int_bt_k_gT, *this->internal_heat_rate, type, ghost_type, -1); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real el_size; Real min_el_size = std::numeric_limits::max(); Real conductivitymax = conductivity(0, 0); // get the biggest parameter from k11 until k33// for (UInt i = 0; i < spatial_dimension; i++) for (UInt j = 0; j < spatial_dimension; j++) conductivitymax = std::max(conductivity(i, j), conductivitymax); - for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { + for (auto & type : + mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array coord(0, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, _not_ghost); auto el_coord = coord.begin(spatial_dimension, nb_nodes_per_element); UInt nb_element = mesh.getNbElement(type); for (UInt el = 0; el < nb_element; ++el, ++el_coord) { el_size = getFEEngine().getElementInradius(*el_coord, type); min_el_size = std::min(min_el_size, el_size); } - + AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size << " and the max conductivity is : " << conductivitymax); } - Real min_dt = - 2. * min_el_size * min_el_size / 4. * density * capacity / conductivitymax; - + Real min_dt = 2. * min_el_size * min_el_size / 4. * density * capacity / + conductivitymax; + mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("heat_transfer").setTimeStep(time_step); #endif } - /* -------------------------------------------------------------------------- */ void HeatTransferModel::readMaterials() { auto sect = this->getParserSection(); if (not std::get<1>(sect)) { const auto & section = std::get<0>(sect); this->parseSection(section); } conductivity_on_qpoints.set(conductivity); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); readMaterials(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacity() { AKANTU_DEBUG_IN(); auto ghost_type = _not_ghost; this->getDOFManager().clearMatrix("M"); auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor rho_functor(*this); - for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { + for (auto && type : + mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldMatrix(rho_functor, "M", "temperature", this->getDOFManager(), type, ghost_type); } need_to_reassemble_capacity = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeRho(Array & rho, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); FEEngine & fem = this->getFEEngine(); UInt nb_element = mesh.getNbElement(type, ghost_type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); rho.resize(nb_element * nb_quadrature_points); rho.set(this->capacity); // Real * rho_1_val = rho.storage(); // /// compute @f$ rho @f$ for each nodes of each element // for (UInt el = 0; el < nb_element; ++el) { // for (UInt n = 0; n < nb_quadrature_points; ++n) { // *rho_1_val++ = this->capacity; // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::computeThermalEnergyByNode() { AKANTU_DEBUG_IN(); Real ethermal = 0.; for (auto && pair : enumerate(make_view( *internal_heat_rate, internal_heat_rate->getNbComponent()))) { auto n = std::get<0>(pair); auto & heat_rate = std::get<1>(pair); Real heat = 0.; bool is_local_node = mesh.isLocalOrMasterNode(n); - //bool is_not_pbc_slave_node = !isPBCSlaveNode(n); - bool count_node = is_local_node;// && is_not_pbc_slave_node; + // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); + bool count_node = is_local_node; // && is_not_pbc_slave_node; for (UInt i = 0; i < heat_rate.size(); ++i) { if (count_node) heat += heat_rate[i] * time_step; } ethermal += heat; } mesh.getCommunicator().allReduce(ethermal, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ethermal; } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::getThermalEnergy( iterator Eth, Array::const_iterator T_it, Array::const_iterator T_end) const { for (; T_it != T_end; ++T_it, ++Eth) { *Eth = capacity * density * *T_it; } } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Vector Eth_on_quarature_points(nb_quadrature_points); auto T_it = this->temperature_on_qpoints(type).begin(); T_it += index * nb_quadrature_points; auto T_end = T_it + nb_quadrature_points; getThermalEnergy(Eth_on_quarature_points.storage(), T_it, T_end); return getFEEngine().integrate(Eth_on_quarature_points, type, index); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy() { Real Eth = 0; auto & fem = getFEEngine(); - for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { + for (auto && type : + mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_element = mesh.getNbElement(type, _not_ghost); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, _not_ghost); Array Eth_per_quad(nb_element * nb_quadrature_points, 1); auto & temperature_interpolated = temperature_on_qpoints(type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type); auto T_it = temperature_interpolated.begin(); auto T_end = temperature_interpolated.end(); getThermalEnergy(Eth_per_quad.begin(), T_it, T_end); Eth += fem.integrate(Eth_per_quad, type); } return Eth; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id) { AKANTU_DEBUG_IN(); Real energy = 0; if (id == "thermal") energy = getThermalEnergy(); // reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id, const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real energy = 0.; if (id == "thermal") energy = getThermalEnergy(type, index); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER -dumper::Field * HeatTransferModel::createNodalFieldBool( +std::shared_ptr HeatTransferModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; - dumper::Field * field = nullptr; + std::shared_ptr field; field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ -dumper::Field * HeatTransferModel::createNodalFieldReal( +std::shared_ptr HeatTransferModel::createNodalFieldReal( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { - if (field_name == "capacity_lumped"){ - AKANTU_EXCEPTION("Capacity lumped is a nodal field now stored in the DOF manager." - "Therefore it cannot be used by a dumper anymore"); + if (field_name == "capacity_lumped") { + AKANTU_EXCEPTION( + "Capacity lumped is a nodal field now stored in the DOF manager." + "Therefore it cannot be used by a dumper anymore"); } std::map *> real_nodal_fields; real_nodal_fields["temperature"] = temperature; real_nodal_fields["temperature_rate"] = temperature_rate; real_nodal_fields["external_heat_rate"] = external_heat_rate; real_nodal_fields["internal_heat_rate"] = internal_heat_rate; real_nodal_fields["increment"] = increment; - dumper::Field * field = + std::shared_ptr field = mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ -dumper::Field * HeatTransferModel::createElementalField( +std::shared_ptr HeatTransferModel::createElementalField( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const UInt & spatial_dimension, const ElementKind & element_kind) { - dumper::Field * field = nullptr; + std::shared_ptr field; if (field_name == "partitions") field = mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); else if (field_name == "temperature_gradient") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(temperature_gradient, element_kind); field = mesh.createElementalField( temperature_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "conductivity") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(conductivity_on_qpoints, element_kind); field = mesh.createElementalField( conductivity_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ -dumper::Field * HeatTransferModel::createElementalField( +std::shared_ptr HeatTransferModel::createElementalField( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const ElementKind & element_kind) { return nullptr; } /* -------------------------------------------------------------------------- */ -dumper::Field * HeatTransferModel::createNodalFieldBool( +std::shared_ptr HeatTransferModel::createNodalFieldBool( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } /* -------------------------------------------------------------------------- */ -dumper::Field * HeatTransferModel::createNodalFieldReal( +std::shared_ptr HeatTransferModel::createNodalFieldReal( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void HeatTransferModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbData(const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = indexes.size(); switch (tag) { case _gst_htm_temperature: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case _gst_htm_temperature: { buffer << (*temperature)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case _gst_htm_temperature: { buffer >> (*temperature)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt HeatTransferModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; Array::const_iterator it = elements.begin(); Array::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case _gst_htm_temperature: { size += nb_nodes_per_element * sizeof(Real); // temperature break; } case _gst_htm_gradient_temperature: { // temperature gradient size += getNbIntegrationPoints(elements) * spatial_dimension * sizeof(Real); size += nb_nodes_per_element * sizeof(Real); // nodal temperatures break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { switch (tag) { case _gst_htm_temperature: { packNodalDataHelper(*temperature, buffer, elements, mesh); break; } case _gst_htm_gradient_temperature: { packElementalDataHelper(temperature_gradient, buffer, elements, true, getFEEngine()); packNodalDataHelper(*temperature, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { switch (tag) { case _gst_htm_temperature: { unpackNodalDataHelper(*temperature, buffer, elements, mesh); break; } case _gst_htm_gradient_temperature: { unpackElementalDataHelper(temperature_gradient, buffer, elements, true, getFEEngine()); unpackNodalDataHelper(*temperature, buffer, elements, mesh); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ -} // akantu +} // namespace akantu diff --git a/src/model/heat_transfer/heat_transfer_model.hh b/src/model/heat_transfer/heat_transfer_model.hh index f96274ada..623a83abd 100644 --- a/src/model/heat_transfer/heat_transfer_model.hh +++ b/src/model/heat_transfer/heat_transfer_model.hh @@ -1,340 +1,343 @@ /** * @file heat_transfer_model.hh * * @author Guillaume Anciaux * @author Lucas Frerot * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Mon Feb 05 2018 * * @brief Model of Heat Transfer * * @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 "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_HEAT_TRANSFER_MODEL_HH__ #define __AKANTU_HEAT_TRANSFER_MODEL_HH__ namespace akantu { template class IntegratorGauss; template class ShapeLagrange; -} +} // namespace akantu namespace akantu { class HeatTransferModel : public Model, public DataAccessor, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using FEEngineType = FEEngineTemplate; HeatTransferModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "heat_transfer_model", const MemoryID & memory_id = 0); virtual ~HeatTransferModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// generic function to initialize everything ready for explicit dynamics void initFullImpl(const ModelOptions & options) override; /// read one material file to instantiate all the materials void readMaterials(); /// allocate all vectors void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize the model void initModel() override; void predictor() override; /// compute the heat flux void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID &) override; /// callback to assemble a Matrix void assembleMatrix(const ID &) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID &) override; std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /* ------------------------------------------------------------------------ */ /* Methods for explicit */ /* ------------------------------------------------------------------------ */ public: /// compute and get the stable time step Real getStableTimeStep(); /// set the stable timestep - void setTimeStep(Real time_step, const ID & solver_id="") override; - -// temporary protection to prevent bad usage: should check for bug + void setTimeStep(Real time_step, const ID & solver_id = "") override; + + // temporary protection to prevent bad usage: should check for bug protected: - /// compute the internal heat flux \todo Need code review: currently not public method + /// compute the internal heat flux \todo Need code review: currently not + /// public method void assembleInternalHeatRate(); public: /// calculate the lumped capacity vector for heat transfer problem void assembleCapacityLumped(); /* ------------------------------------------------------------------------ */ /* Methods for static */ /* ------------------------------------------------------------------------ */ public: /// assemble the conductivity matrix void assembleConductivityMatrix(); /// assemble the conductivity matrix void assembleCapacity(); /// compute the capacity on quadrature points void computeRho(Array & rho, ElementType type, GhostType ghost_type); private: /// calculate the lumped capacity vector for heat transfer problem (w /// ghost type) void assembleCapacityLumped(const GhostType & ghost_type); /// assemble the conductivity matrix (w/ ghost type) template void assembleConductivityMatrix(const GhostType & ghost_type); /// compute the conductivity tensor for each quadrature point in an array void computeConductivityOnQuadPoints(const GhostType & ghost_type); /// compute vector k \grad T for each quadrature point void computeKgradT(const GhostType & ghost_type); /// compute the thermal energy Real computeThermalEnergyByNode(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; inline UInt getNbData(const Array & indexes, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: - dumper::Field * createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - dumper::Field * createNodalFieldBool(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - dumper::Field * createElementalField(const std::string & field_name, - const std::string & group_name, - bool padding_flag, - const UInt & spatial_dimension, - const ElementKind & kind) override; + std::shared_ptr + createNodalFieldReal(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; + + std::shared_ptr + createNodalFieldBool(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; + + std::shared_ptr + createElementalField(const std::string & field_name, + const std::string & group_name, bool padding_flag, + const UInt & spatial_dimension, + const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Density, density, Real); AKANTU_GET_MACRO(Capacity, capacity, Real); /// get the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// get the assembled heat flux AKANTU_GET_MACRO(InternalHeatRate, *internal_heat_rate, Array &); /// get the boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the external heat rate vector AKANTU_GET_MACRO(ExternalHeatRate, *external_heat_rate, Array &); /// get the temperature gradient AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureGradient, temperature_gradient, Real); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConductivityOnQpoints, conductivity_on_qpoints, Real); /// get the conductivity on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureOnQpoints, temperature_on_qpoints, Real); /// internal variables AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(KgradT, k_gradt_on_qpoints, Real); /// get the temperature AKANTU_GET_MACRO(Temperature, *temperature, Array &); /// get the temperature derivative AKANTU_GET_MACRO(TemperatureRate, *temperature_rate, Array &); /// get the energy denominated by thermal Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); /// get the energy denominated by thermal Real getEnergy(const std::string & energy_id); /// get the thermal energy for a given element Real getThermalEnergy(const ElementType & type, UInt index); /// get the thermal energy for a given element Real getThermalEnergy(); protected: /* ------------------------------------------------------------------------ */ FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ----------------------------------------------------------------------- */ template void getThermalEnergy(iterator Eth, Array::const_iterator T_it, Array::const_iterator T_end) const; template void allocNodalField(Array *& array, const ID & name); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// number of iterations UInt n_iter; /// time step Real time_step; /// temperatures array Array * temperature{nullptr}; /// temperatures derivatives array Array * temperature_rate{nullptr}; /// increment array (@f$\delta \dot T@f$ or @f$\delta T@f$) Array * increment{nullptr}; /// the density Real density; /// the speed of the changing temperature ElementTypeMapArray temperature_gradient; /// temperature field on quadrature points ElementTypeMapArray temperature_on_qpoints; /// conductivity tensor on quadrature points ElementTypeMapArray conductivity_on_qpoints; /// vector k \grad T on quad points ElementTypeMapArray k_gradt_on_qpoints; /// external flux vector Array * external_heat_rate{nullptr}; /// residuals array Array * internal_heat_rate{nullptr}; /// boundary vector Array * blocked_dofs{nullptr}; // realtime Real time; /// capacity Real capacity; // conductivity matrix Matrix conductivity; // linear variation of the conductivity (for temperature dependent // conductivity) Real conductivity_variation; // reference temperature for the interpretation of temperature variation Real T_ref; // the biggest parameter of conductivity matrix Real conductivitymax; bool need_to_reassemble_capacity{true}; bool need_to_reassemble_capacity_lumped{true}; UInt temperature_release{0}; UInt conductivity_matrix_release{0}; std::unordered_map initial_conductivity{{_not_ghost, true}, {_ghost, true}}; std::unordered_map conductivity_release{{_not_ghost, 0}, {_ghost, 0}}; }; -} // akantu +} // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_model_inline_impl.cc" #endif /* __AKANTU_HEAT_TRANSFER_MODEL_HH__ */ diff --git a/src/model/model.cc b/src/model/model.cc index acd295fd7..085695e05 100644 --- a/src/model/model.cc +++ b/src/model/model.cc @@ -1,371 +1,372 @@ /** * @file model.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Mon Oct 03 2011 * @date last modification: Tue Feb 20 2018 * * @brief implementation of model common parts * * @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 "model.hh" #include "communicator.hh" #include "data_accessor.hh" #include "element_group.hh" #include "element_synchronizer.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Model::Model(Mesh & mesh, const ModelType & type, UInt dim, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), ModelSolver(mesh, type, id, memory_id), mesh(mesh), spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension() : dim), parser(getStaticParser()) { AKANTU_DEBUG_IN(); this->mesh.registerEventHandler(*this, _ehp_model); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Model::~Model() = default; /* -------------------------------------------------------------------------- */ // void Model::setParser(Parser & parser) { this->parser = &parser; } /* -------------------------------------------------------------------------- */ void Model::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); method = options.analysis_method; if (!this->hasDefaultSolver()) { this->initNewSolver(this->method); } initModel(); initFEEngineBoundary(); - //if(mesh.isPeriodic()) this->initPBC(); + // if(mesh.isPeriodic()) this->initPBC(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Model::initNewSolver(const AnalysisMethod & method) { ID solver_name; TimeStepSolverType tss_type; std::tie(solver_name, tss_type) = this->getDefaultSolverID(method); if (not this->hasSolver(solver_name)) { ModelSolverOptions options = this->getDefaultSolverOptions(tss_type); this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type); for (auto && is_type : options.integration_scheme_type) { if (!this->hasIntegrationScheme(solver_name, is_type.first)) { this->setIntegrationScheme(solver_name, is_type.first, is_type.second, options.solution_type[is_type.first]); } } } this->method = method; this->setDefaultSolver(solver_name); } /* -------------------------------------------------------------------------- */ // void Model::initPBC() { // auto it = pbc_pair.begin(); // auto end = pbc_pair.end(); // is_pbc_slave_node.resize(mesh.getNbNodes()); //#ifndef AKANTU_NDEBUG // auto coord_it = mesh.getNodes().begin(this->spatial_dimension); //#endif // while (it != end) { // UInt i1 = (*it).first; // is_pbc_slave_node(i1) = true; // #ifndef AKANTU_NDEBUG // UInt i2 = (*it).second; // UInt slave = mesh.isDistributed() ? mesh.getGlobalNodesIds()(i1) : i1; // UInt master = mesh.isDistributed() ? mesh.getGlobalNodesIds()(i2) : i2; -// AKANTU_DEBUG_INFO("pairing " << slave << " (" << Vector(coord_it[i1]) +// AKANTU_DEBUG_INFO("pairing " << slave << " (" << +// Vector(coord_it[i1]) // << ") with " << master << " (" // << Vector(coord_it[i2]) << ")"); // #endif // ++it; // } // } /* -------------------------------------------------------------------------- */ void Model::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup(const std::string & group_name, const std::string & dumper_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void Model::dumpGroup() { auto bit = mesh.element_group_begin(); auto bend = mesh.element_group_end(); for (; bit != bend; ++bit) { bit->second->dump(); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory) { auto bit = mesh.element_group_begin(); auto bend = mesh.element_group_end(); for (; bit != bend; ++bit) { bit->second->setDirectory(directory); } } /* -------------------------------------------------------------------------- */ void Model::setGroupDirectory(const std::string & directory, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Model::setGroupBaseName(const std::string & basename, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.setBaseName(basename); } /* -------------------------------------------------------------------------- */ DumperIOHelper & Model::getGroupDumper(const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); return group.getDumper(); } /* -------------------------------------------------------------------------- */ // DUMPER stuff /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & field_id, - dumper::Field * field, + std::shared_ptr field, DumperIOHelper & dumper) { #ifdef AKANTU_USE_IOHELPER dumper.registerField(field_id, field); #endif } /* -------------------------------------------------------------------------- */ void Model::addDumpField(const std::string & field_id) { this->addDumpFieldToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldVector(const std::string & field_id) { this->addDumpFieldVectorToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldTensor(const std::string & field_id) { this->addDumpFieldTensorToDumper(mesh.getDefaultDumperName(), field_id); } /* -------------------------------------------------------------------------- */ void Model::setBaseName(const std::string & field_id) { mesh.setBaseName(field_id); } /* -------------------------------------------------------------------------- */ void Model::setBaseNameToDumper(const std::string & dumper_name, const std::string & basename) { mesh.setBaseNameToDumper(dumper_name, basename); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular, false); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id, group_name, _ek_regular, false); } /* -------------------------------------------------------------------------- */ void Model::removeDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void Model::removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.removeDumpFieldFromDumper(dumper_name, field_id); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name, _ek_regular, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular, true); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name, this->spatial_dimension, element_kind, padding_flag); } /* -------------------------------------------------------------------------- */ void Model::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, const ElementKind & element_kind, bool padding_flag) { #ifdef AKANTU_USE_IOHELPER - dumper::Field * field = nullptr; + std::shared_ptr field; if (!field) field = this->createNodalFieldReal(field_id, group_name, padding_flag); if (!field) field = this->createNodalFieldUInt(field_id, group_name, padding_flag); if (!field) field = this->createNodalFieldBool(field_id, group_name, padding_flag); if (!field) field = this->createElementalField(field_id, group_name, padding_flag, spatial_dimension, element_kind); if (!field) field = this->mesh.createFieldFromAttachedData(field_id, group_name, element_kind); if (!field) field = this->mesh.createFieldFromAttachedData(field_id, group_name, element_kind); #ifndef AKANTU_NDEBUG if (!field) { AKANTU_DEBUG_WARNING("No field could be found based on name: " << field_id); } #endif if (field) { DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name, group_name); this->addDumpGroupFieldToDumper(field_id, field, dumper); } #endif } /* -------------------------------------------------------------------------- */ void Model::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void Model::setDirectory(const std::string & directory) { mesh.setDirectory(directory); } /* -------------------------------------------------------------------------- */ void Model::setDirectoryToDumper(const std::string & dumper_name, const std::string & directory) { mesh.setDirectoryToDumper(dumper_name, directory); } /* -------------------------------------------------------------------------- */ void Model::setTextModeToDumper() { mesh.setTextModeToDumper(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/model.hh b/src/model/model.hh index 01eb8e162..3f81f4c14 100644 --- a/src/model/model.hh +++ b/src/model/model.hh @@ -1,344 +1,344 @@ /** * @file model.hh * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief Interface of a model * * @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_common.hh" #include "aka_memory.hh" #include "aka_named_argument.hh" #include "fe_engine.hh" #include "mesh.hh" #include "model_options.hh" #include "model_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MODEL_HH__ #define __AKANTU_MODEL_HH__ namespace akantu { class SynchronizerRegistry; class Parser; class DumperIOHelper; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { class Model : public Memory, public ModelSolver, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Model(Mesh & mesh, const ModelType & type, UInt spatial_dimension = _all_dimensions, const ID & id = "model", const MemoryID & memory_id = 0); ~Model() override; using FEEngineMap = std::map>; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: virtual void initFullImpl(const ModelOptions & options); public: template std::enable_if_t::value> initFull(pack &&... _pack) { switch (this->model_type) { #ifdef AKANTU_SOLID_MECHANICS case ModelType::_solid_mechanics_model: this->initFullImpl(SolidMechanicsModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_COHESIVE_ELEMENT case ModelType::_solid_mechanics_model_cohesive: this->initFullImpl(SolidMechanicsModelCohesiveOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_HEAT_TRANSFER case ModelType::_heat_transfer_model: this->initFullImpl(HeatTransferModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif #ifdef AKANTU_EMBEDDED case ModelType::_embedded_model: this->initFullImpl(EmbeddedInterfaceModelOptions{ use_named_args, std::forward(_pack)...}); break; #endif default: this->initFullImpl(ModelOptions{use_named_args, std::forward(_pack)...}); } } template std::enable_if_t::value> initFull(pack &&... _pack) { this->initFullImpl(std::forward(_pack)...); } /// initialize a new solver if needed void initNewSolver(const AnalysisMethod & method); protected: /// get some default values for derived classes virtual std::tuple getDefaultSolverID(const AnalysisMethod & method) = 0; virtual void initModel() = 0; virtual void initFEEngineBoundary(); /// function to print the containt of the class void printself(std::ostream &, int = 0) const override{}; public: /* ------------------------------------------------------------------------ */ /* Access to the dumpable interface of the boundaries */ /* ------------------------------------------------------------------------ */ /// Dump the data for a given group void dumpGroup(const std::string & group_name); void dumpGroup(const std::string & group_name, const std::string & dumper_name); /// Dump the data for all boundaries void dumpGroup(); /// Set the directory for a given group void setGroupDirectory(const std::string & directory, const std::string & group_name); /// Set the directory for all boundaries void setGroupDirectory(const std::string & directory); /// Set the base name for a given group void setGroupBaseName(const std::string & basename, const std::string & group_name); /// Get the internal dumper of a given group DumperIOHelper & getGroupDumper(const std::string & group_name); /* ------------------------------------------------------------------------ */ /* Function for non local capabilities */ /* ------------------------------------------------------------------------ */ virtual void updateDataForNonLocalCriterion(__attribute__((unused)) ElementTypeMapReal & criterion) { AKANTU_TO_IMPLEMENT(); } protected: template void allocNodalField(Array *& array, UInt nb_component, const ID & name); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get id of model AKANTU_GET_MACRO(ID, id, const ID) /// get the number of surfaces AKANTU_GET_MACRO(Mesh, mesh, Mesh &) /// synchronize the boundary in case of parallel run virtual void synchronizeBoundaries(){}; /// return the fem object associated with a provided name inline FEEngine & getFEEngine(const ID & name = "") const; /// return the fem boundary object associated with a provided name virtual FEEngine & getFEEngineBoundary(const ID & name = ""); /// register a fem object associated with name template inline void registerFEEngineObject(const std::string & name, Mesh & mesh, UInt spatial_dimension); /// unregister a fem object associated with name inline void unRegisterFEEngineObject(const std::string & name); /// return the synchronizer registry SynchronizerRegistry & getSynchronizerRegistry(); /// return the fem object associated with a provided name template inline FEEngineClass & getFEEngineClass(std::string name = "") const; /// return the fem boundary object associated with a provided name template inline FEEngineClass & getFEEngineClassBoundary(std::string name = ""); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /* ------------------------------------------------------------------------ */ /* Pack and unpack helper functions */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbIntegrationPoints(const Array & elements, const ID & fem_id = ID()) const; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ void setTextModeToDumper(); virtual void addDumpGroupFieldToDumper(const std::string & field_id, - dumper::Field * field, + std::shared_ptr field, DumperIOHelper & dumper); virtual void addDumpField(const std::string & field_id); virtual void addDumpFieldVector(const std::string & field_id); virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpFieldTensor(const std::string & field_id); virtual void setBaseName(const std::string & basename); virtual void setBaseNameToDumper(const std::string & dumper_name, const std::string & basename); virtual void addDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, UInt spatial_dimension, const ElementKind & element_kind, bool padding_flag); virtual void removeDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); - virtual dumper::Field * + virtual std::shared_ptr createNodalFieldReal(__attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } - virtual dumper::Field * + virtual std::shared_ptr createNodalFieldUInt(__attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } - virtual dumper::Field * + virtual std::shared_ptr createNodalFieldBool(__attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } - virtual dumper::Field * + virtual std::shared_ptr createElementalField(__attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag, __attribute__((unused)) const UInt & spatial_dimension, __attribute__((unused)) const ElementKind & kind) { return nullptr; } void setDirectory(const std::string & directory); void setDirectoryToDumper(const std::string & dumper_name, const std::string & directory); virtual void dump(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: friend std::ostream & operator<<(std::ostream &, const Model &); /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// Mesh Mesh & mesh; /// Spatial dimension of the problem UInt spatial_dimension; /// the main fem object present in all models FEEngineMap fems; /// the fem object present in all models for boundaries FEEngineMap fems_boundary; /// default fem object std::string default_fem; /// parser to the pointer to use Parser & parser; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Model & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "model_inline_impl.cc" #endif /* __AKANTU_MODEL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 486414299..3f95aa752 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,559 +1,561 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * @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 "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager_callback.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel( Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep() override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(const GhostType & ghost_type) override; void insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, const ElementKind & element_kind); //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, const ElementKind & kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); - dumper::Field * createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - dumper::Field * createNodalFieldBool(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - dumper::Field * createElementalField(const std::string & field_name, - const std::string & group_name, - bool padding_flag, - const UInt & spatial_dimension, - const ElementKind & kind) override; + std::shared_ptr + createNodalFieldReal(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; + + std::shared_ptr + createNodalFieldBool(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; + + std::shared_ptr + createElementalField(const std::string & field_name, + const std::string & group_name, bool padding_flag, + const UInt & spatial_dimension, + const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition const Array & getCurrentPosition(); /// get the SolidMechanicsModel::increment vector \warn only consistent if AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the SolidMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the SolidMechanicsModel::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the SolidMechanicsModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector, MaterialSelector &); AKANTU_SET_MACRO(MaterialSelector, material_selector, std::shared_ptr); /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; UInt displacement_release{0}; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// lumped mass array Array * mass{nullptr}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// velocities array Array * velocity{nullptr}; /// accelerations array Array * acceleration{nullptr}; /// accelerations array // Array * increment_acceleration; /// external forces array Array * external_force{nullptr}; /// internal forces array Array * internal_force{nullptr}; /// array specifing if a degree of freedom is blocked or not Array * blocked_dofs{nullptr}; /// array of current position used during update residual Array * current_position{nullptr}; UInt current_position_release{0}; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; /// list of used materials std::vector> materials; /// mapping between material name and material internal id std::map materials_names_to_id; /// class defining of to choose a material std::shared_ptr material_selector; /// tells if the material are instantiated bool are_materials_instantiated; using flatten_internal_map = std::map, ElementTypeMapArray *>; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.cc" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc index 4c260b2df..bead9afbe 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_embedded_interface/embedded_interface_model.cc @@ -1,173 +1,173 @@ /** * @file embedded_interface_model.cc * * @author Lucas Frerot * * @date creation: Fri Mar 13 2015 * @date last modification: Wed Feb 14 2018 * * @brief Model of Solid Mechanics with embedded interfaces * * @section LICENSE * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "embedded_interface_model.hh" #include "integrator_gauss.hh" #include "material_elastic.hh" #include "material_reinforcement.hh" #include "mesh_iterators.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumpable_inline_impl.hh" #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ EmbeddedInterfaceModel::EmbeddedInterfaceModel(Mesh & mesh, Mesh & primitive_mesh, UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, spatial_dimension, id, memory_id), intersector(mesh, primitive_mesh), interface_mesh(nullptr), primitive_mesh(primitive_mesh), interface_material_selector(nullptr) { this->model_type = ModelType::_embedded_model; // This pointer should be deleted by ~SolidMechanicsModel() auto mat_sel_pointer = std::make_shared>("physical_names", *this); this->setMaterialSelector(mat_sel_pointer); interface_mesh = &(intersector.getInterfaceMesh()); // Create 1D FEEngine on the interface mesh registerFEEngineObject("EmbeddedInterfaceFEEngine", *interface_mesh, 1); // Registering allocator for material reinforcement MaterialFactory::getInstance().registerAllocator( "reinforcement", [&](UInt dim, const ID & constitutive, SolidMechanicsModel &, const ID & id) -> std::unique_ptr { if (constitutive == "elastic") { using mat = MaterialElastic<1>; switch (dim) { case 2: return std::make_unique>(*this, id); case 3: return std::make_unique>(*this, id); default: AKANTU_EXCEPTION("Dimension 1 is invalid for reinforcements"); } } else { AKANTU_EXCEPTION("Reinforcement type" << constitutive << " is not recognized"); } }); } /* -------------------------------------------------------------------------- */ EmbeddedInterfaceModel::~EmbeddedInterfaceModel() { delete interface_material_selector; } /* -------------------------------------------------------------------------- */ void EmbeddedInterfaceModel::initFullImpl(const ModelOptions & options) { const auto & eim_options = dynamic_cast(options); // Do no initialize interface_mesh if told so if (eim_options.has_intersections) intersector.constructData(); SolidMechanicsModel::initFullImpl(options); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("reinforcement", id); this->mesh.addDumpMeshToDumper("reinforcement", *interface_mesh, 1, _not_ghost, _ek_regular); #endif } void EmbeddedInterfaceModel::initModel() { // Initialize interface FEEngine SolidMechanicsModel::initModel(); FEEngine & engine = getFEEngine("EmbeddedInterfaceFEEngine"); engine.initShapeFunctions(_not_ghost); engine.initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void EmbeddedInterfaceModel::assignMaterialToElements( const ElementTypeMapArray * filter) { delete interface_material_selector; interface_material_selector = new InterfaceMeshDataMaterialSelector("physical_names", *this); for_each_element(getInterfaceMesh(), [&](auto && element) { auto mat_index = (*interface_material_selector)(element); // material_index(element) = mat_index; materials[mat_index]->addElement(element); // this->material_local_numbering(element) = index; }, _element_filter = filter, _spatial_dimension = 1); SolidMechanicsModel::assignMaterialToElements(filter); } /* -------------------------------------------------------------------------- */ void EmbeddedInterfaceModel::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { #ifdef AKANTU_USE_IOHELPER - dumper::Field * field = NULL; + std::shared_ptr field; // If dumper is reinforcement, create a 1D elemental field if (dumper_name == "reinforcement") field = this->createElementalField(field_id, group_name, padding_flag, 1, element_kind); else { try { SolidMechanicsModel::addDumpGroupFieldToDumper( dumper_name, field_id, group_name, element_kind, padding_flag); } catch (...) { } } if (field) { DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name, group_name); Model::addDumpGroupFieldToDumper(field_id, field, dumper); } #endif } } // akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_io.cc b/src/model/solid_mechanics/solid_mechanics_model_io.cc index ba29a79ca..b371b0373 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_io.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_io.cc @@ -1,326 +1,325 @@ /** * @file solid_mechanics_model_io.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Sun Jul 09 2017 * @date last modification: Sun Dec 03 2017 * * @brief Dumpable part of the SolidMechnicsModel * * @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 "solid_mechanics_model.hh" #include "group_manager_inline_impl.cc" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_field.hh" #include "dumper_homogenizing_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper.hh" #include "dumper_material_padders.hh" #include "dumper_paraview.hh" #endif namespace akantu { /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::isInternal(const std::string & field_name, const ElementKind & element_kind) { /// check if at least one material contains field_id as an internal for (auto & material : materials) { bool is_internal = material->isInternal(field_name, element_kind); if (is_internal) return true; } return false; } /* -------------------------------------------------------------------------- */ ElementTypeMap SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name, const ElementKind & element_kind) { if (!(this->isInternal(field_name, element_kind))) AKANTU_EXCEPTION("unknown internal " << field_name); for (auto & material : materials) { if (material->isInternal(field_name, element_kind)) return material->getInternalDataPerElem(field_name, element_kind); } return ElementTypeMap(); } /* -------------------------------------------------------------------------- */ ElementTypeMapArray & SolidMechanicsModel::flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type) { std::pair key(field_name, kind); if (this->registered_internals.count(key) == 0) { this->registered_internals[key] = new ElementTypeMapArray(field_name, this->id, this->memory_id); } ElementTypeMapArray * internal_flat = this->registered_internals[key]; for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, kind)) { if (internal_flat->exists(type, ghost_type)) { auto & internal = (*internal_flat)(type, ghost_type); // internal.clear(); internal.resize(0); } } for (auto & material : materials) { if (material->isInternal(field_name, kind)) material->flattenInternal(field_name, *internal_flat, ghost_type, kind); } return *internal_flat; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::flattenAllRegisteredInternals( const ElementKind & kind) { ElementKind _kind; ID _id; for (auto & internal : this->registered_internals) { std::tie(_id, _kind) = internal.first; if (kind == _kind) this->flattenInternal(_id, kind); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onDump() { this->flattenAllRegisteredInternals(_ek_regular); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER -dumper::Field * SolidMechanicsModel::createElementalField( +std::shared_ptr SolidMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { - dumper::Field * field = nullptr; + std::shared_ptr field; if (field_name == "partitions") field = mesh.createElementalField( mesh.getConnectivities(), group_name, spatial_dimension, kind); else if (field_name == "material_index") field = mesh.createElementalField( material_index, group_name, spatial_dimension, kind); else { // this copy of field_name is used to compute derivated data such as // strain and von mises stress that are based on grad_u and stress std::string field_name_copy(field_name); if (field_name == "strain" || field_name == "Green strain" || field_name == "principal strain" || field_name == "principal Green strain") field_name_copy = "grad_u"; else if (field_name == "Von Mises stress") field_name_copy = "stress"; bool is_internal = this->isInternal(field_name_copy, kind); if (is_internal) { ElementTypeMap nb_data_per_elem = this->getInternalDataPerElem(field_name_copy, kind); ElementTypeMapArray & internal_flat = this->flattenInternal(field_name_copy, kind); field = mesh.createElementalField( internal_flat, group_name, spatial_dimension, kind, nb_data_per_elem); if (field_name == "strain") { auto * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "Von Mises stress") { auto * foo = new dumper::ComputeVonMisesStress(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "Green strain") { auto * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "principal strain") { auto * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "principal Green strain") { auto * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } // treat the paddings if (padding_flag) { if (field_name == "stress") { if (spatial_dimension == 2) { auto * foo = new dumper::StressPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } } else if (field_name == "strain" || field_name == "Green strain") { if (spatial_dimension == 2) { auto * foo = new dumper::StrainPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } } } // homogenize the field - dumper::ComputeFunctorInterface * foo = - dumper::HomogenizerProxy::createHomogenizer(*field); + auto foo = dumper::HomogenizerProxy::createHomogenizer(*field); field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } } return field; } /* -------------------------------------------------------------------------- */ -dumper::Field * +std::shared_ptr SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["displacement"] = this->displacement; real_nodal_fields["mass"] = this->mass; real_nodal_fields["velocity"] = this->velocity; real_nodal_fields["acceleration"] = this->acceleration; real_nodal_fields["external_force"] = this->external_force; real_nodal_fields["internal_force"] = this->internal_force; real_nodal_fields["increment"] = this->displacement_increment; if (field_name == "force") { AKANTU_EXCEPTION("The 'force' field has been renamed in 'external_force'"); } else if (field_name == "residual") { AKANTU_EXCEPTION( "The 'residual' field has been replaced by 'internal_force'"); } - dumper::Field * field = nullptr; + std::shared_ptr field; if (padding_flag) field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); else field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ -dumper::Field * SolidMechanicsModel::createNodalFieldBool( +std::shared_ptr SolidMechanicsModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; - dumper::Field * field = nullptr; + std::shared_ptr field; field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ -dumper::Field * SolidMechanicsModel::createElementalField(const std::string &, - const std::string &, - bool, const UInt &, - const ElementKind &) { +std::shared_ptr +SolidMechanicsModel::createElementalField(const std::string &, + const std::string &, bool, + const UInt &, const ElementKind &) { return nullptr; } /* -------------------------------------------------------------------------- */ -dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string &, - const std::string &, - bool) { +std::shaed_ptr +SolidMechanicsModel::createNodalFieldReal(const std::string &, + const std::string &, bool) { return nullptr; } /* -------------------------------------------------------------------------- */ -dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string &, - const std::string &, - bool) { +std::shared_ptr +SolidMechanicsModel::createNodalFieldBool(const std::string &, + const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump() { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(time, step); } } // namespace akantu diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index 91d1bc1d8..20a69bbb6 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,457 +1,455 @@ /** * @file structural_mechanics_model.cc * * @author Fabian Barras * @author Lucas Frerot * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Wed Feb 21 2018 * * @brief Model implementation for StucturalMechanics elements * * @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 "structural_mechanics_model.hh" #include "dof_manager.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "shape_structural.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumpable_inline_impl.hh" #include "dumper_elemental_field.hh" #include "dumper_iohelper_paraview.hh" #include "group_manager_inline_impl.cc" #endif /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model_inline_impl.cc" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, ModelType::_structural_mechanics_model, dim, id, memory_id), time_step(NAN), f_m2a(1.0), stress("stress", id, memory_id), element_material("element_material", id, memory_id), set_ID("beam sets", id, memory_id), rotation_matrix("rotation_matices", id, memory_id) { AKANTU_DEBUG_IN(); registerFEEngineObject("StructuralMechanicsFEEngine", mesh, spatial_dimension); if (spatial_dimension == 2) nb_degree_of_freedom = 3; else if (spatial_dimension == 3) nb_degree_of_freedom = 6; else { AKANTU_TO_IMPLEMENT(); } #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("paraview_all", id, true); #endif this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); this->initDOFManager(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) { // <<<< This is the SolidMechanicsModel implementation for future ref >>>> // material_index.initialize(mesh, _element_kind = _ek_not_defined, // _default_value = UInt(-1), _with_nb_element = // true); // material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, // _with_nb_element = true); // Model::initFullImpl(options); // // initialize pbc // if (this->pbc_pair.size() != 0) // this->initPBC(); // // initialize the materials // if (this->parser.getLastParsedFile() != "") { // this->instantiateMaterials(); // } // this->initMaterials(); // this->initBC(*this, *displacement, *displacement_increment, // *external_force); // <<<< END >>>> Model::initFullImpl(options); // Initializing stresses ElementTypeMap stress_components; /// TODO this is ugly af, maybe add a function to FEEngine for (auto & type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { UInt nb_components = 0; // Getting number of components for each element type #define GET_(type) nb_components = ElementClass::getNbStressComponents() AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_); #undef GET_ stress_components(nb_components, type); } - stress.initialize(getFEEngine(), _spatial_dimension = _all_dimensions, - _element_kind = _ek_structural, _all_ghost_types = true, - _nb_component = [&stress_components]( - const ElementType & type, const GhostType &) -> UInt { - return stress_components(type); - }); + stress.initialize( + getFEEngine(), _spatial_dimension = _all_dimensions, + _element_kind = _ek_structural, _all_ghost_types = true, + _nb_component = [&stress_components](const ElementType & type, + const GhostType &) -> UInt { + return stress_components(type); + }); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFEEngineBoundary() { /// TODO: this function should not be reimplemented /// we're just avoiding a call to Model::initFEEngineBoundary() } /* -------------------------------------------------------------------------- */ // void StructuralMechanicsModel::setTimeStep(Real time_step) { // this->time_step = time_step; // #if defined(AKANTU_USE_IOHELPER) // this->mesh.getDumper().setTimeStep(time_step); // #endif // } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType) { AKANTU_DEBUG_IN(); this->allocNodalField(displacement_rotation, nb_degree_of_freedom, "displacement"); this->allocNodalField(external_force, nb_degree_of_freedom, "external_force"); this->allocNodalField(internal_force, nb_degree_of_freedom, "internal_force"); this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs"); auto & dof_manager = this->getDOFManager(); if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *displacement_rotation, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); } if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(velocity, spatial_dimension, "velocity"); this->allocNodalField(acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { for (auto && type : mesh.elementTypes(_element_kind = _ek_structural)) { // computeRotationMatrix(type); element_material.alloc(mesh.getNbElement(type), 1, type); } getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); getDOFManager().getMatrix("K").clear(); for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define ASSEMBLE_STIFFNESS_MATRIX(type) assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define COMPUTE_STRESS_ON_QUAD(type) computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeRotationMatrix(const ElementType & type) { Mesh & mesh = getFEEngine().getMesh(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); if (!rotation_matrix.exists(type)) { rotation_matrix.alloc(nb_element, nb_degree_of_freedom * nb_nodes_per_element * nb_degree_of_freedom * nb_nodes_per_element, type); } else { rotation_matrix(type).resize(nb_element); } rotation_matrix(type).clear(); Array rotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom); rotations.clear(); #define COMPUTE_ROTATION_MATRIX(type) computeRotationMatrix(rotations); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX); #undef COMPUTE_ROTATION_MATRIX auto R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); auto T_it = rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element, nb_degree_of_freedom * nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) { auto & T = *T_it; auto & R = *R_it; for (UInt k = 0; k < nb_nodes_per_element; ++k) { for (UInt i = 0; i < nb_degree_of_freedom; ++i) for (UInt j = 0; j < nb_degree_of_freedom; ++j) T(k * nb_degree_of_freedom + i, k * nb_degree_of_freedom + j) = R(i, j); } } } /* -------------------------------------------------------------------------- */ -dumper::Field * StructuralMechanicsModel::createNodalFieldBool( +std::shared_ptr StructuralMechanicsModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; - dumper::Field * field = NULL; - field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); - return field; + return mesh.createNodalField(uint_nodal_fields[field_name], group_name); } /* -------------------------------------------------------------------------- */ -dumper::Field * +std::shared_ptr StructuralMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { UInt n; if (spatial_dimension == 2) { n = 2; } else { n = 3; } if (field_name == "displacement") { return mesh.createStridedNodalField(displacement_rotation, group_name, n, 0, padding_flag); } if (field_name == "rotation") { return mesh.createStridedNodalField(displacement_rotation, group_name, nb_degree_of_freedom - n, n, padding_flag); } if (field_name == "force") { return mesh.createStridedNodalField(external_force, group_name, n, 0, padding_flag); } if (field_name == "momentum") { return mesh.createStridedNodalField( external_force, group_name, nb_degree_of_freedom - n, n, padding_flag); } if (field_name == "internal_force") { return mesh.createStridedNodalField(internal_force, group_name, n, 0, padding_flag); } if (field_name == "internal_momentum") { return mesh.createStridedNodalField( internal_force, group_name, nb_degree_of_freedom - n, n, padding_flag); } return nullptr; } /* -------------------------------------------------------------------------- */ -dumper::Field * StructuralMechanicsModel::createElementalField( +std::shared_ptr StructuralMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool, - const UInt & spatial_dimension, - const ElementKind & kind) { + const UInt & spatial_dimension, const ElementKind & kind) { - dumper::Field * field = NULL; + std::shared_ptr field; if (field_name == "element_index_by_material") field = mesh.createElementalField( field_name, group_name, spatial_dimension, kind); return field; } /* -------------------------------------------------------------------------- */ /* Virtual methods from SolverCallback */ /* -------------------------------------------------------------------------- */ /// get the type of matrix needed MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) { return _symmetric; } /// callback to assemble a Matrix void StructuralMechanicsModel::assembleMatrix(const ID & id) { if (id == "K") assembleStiffnessMatrix(); } /// callback to assemble a lumped Matrix void StructuralMechanicsModel::assembleLumpedMatrix(const ID & /*id*/) {} /// callback to assemble the residual StructuralMechanicsModel::(rhs) void StructuralMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); auto & dof_manager = getDOFManager(); internal_force->clear(); computeStresses(); assembleInternalForce(); dof_manager.assembleToResidual("displacement", *internal_force, -1); dof_manager.assembleToResidual("displacement", *external_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Virtual methods from Model */ /* -------------------------------------------------------------------------- */ /// get some default values for derived classes std::tuple StructuralMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _static: { return std::make_tuple("static", _tsst_static); } case _implicit_dynamic: { return std::make_tuple("implicit", _tsst_dynamic); } default: return std::make_tuple("unknown", _tsst_not_defined); } } /* ------------------------------------------------------------------------ */ ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_static: { options.non_linear_solver_type = _nls_linear; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce() { for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { assembleInternalForce(type, _not_ghost); // assembleInternalForce(type, _ghost); } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce(const ElementType & type, GhostType gt) { auto & fem = getFEEngine(); auto & sigma = stress(type, gt); auto ndof = getNbDegreeOfFreedom(type); auto nb_nodes = mesh.getNbNodesPerElement(type); auto ndof_per_elem = ndof * nb_nodes; Array BtSigma(fem.getNbIntegrationPoints(type) * mesh.getNbElement(type), ndof_per_elem, "BtSigma"); fem.computeBtD(sigma, BtSigma, type, gt); Array intBtSigma(0, ndof_per_elem, "intBtSigma"); fem.integrate(BtSigma, intBtSigma, ndof_per_elem, type, gt); BtSigma.resize(0); getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force, type, gt, 1); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh index d7a169f22..d6a3dbcc7 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.hh +++ b/src/model/structural_mechanics/structural_mechanics_model.hh @@ -1,309 +1,311 @@ /** * @file structural_mechanics_model.hh * * @author Fabian Barras * @author Lucas Frerot * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Tue Feb 20 2018 * * @brief Particular implementation of the structural elements in the * StructuralMechanicsModel * * @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_named_argument.hh" #include "boundary_condition.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ #define __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeStructural; } // namespace akantu namespace akantu { struct StructuralMaterial { Real E{0}; Real A{1}; Real I{0}; Real Iz{0}; Real Iy{0}; Real GJ{0}; Real rho{0}; Real t{0}; Real nu{0}; }; class StructuralMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineType = FEEngineTemplate; StructuralMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "structural_mechanics_model", const MemoryID & memory_id = 0); virtual ~StructuralMechanicsModel(); /// Init full model void initFullImpl(const ModelOptions & options) override; /// Init boundary FEEngine void initFEEngineBoundary() override; /* ------------------------------------------------------------------------ */ /* Virtual methods from SolverCallback */ /* ------------------------------------------------------------------------ */ /// get the type of matrix needed MatrixType getMatrixType(const ID &) override; /// callback to assemble a Matrix void assembleMatrix(const ID &) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID &) override; /// callback to assemble the residual (rhs) void assembleResidual() override; /* ------------------------------------------------------------------------ */ /* Virtual methods from Model */ /* ------------------------------------------------------------------------ */ protected: /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; UInt getNbDegreeOfFreedom(const ElementType & type) const; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize the model void initModel() override; /// compute the stresses per elements void computeStresses(); /// compute the nodal forces void assembleInternalForce(); /// compute the nodal forces for an element type void assembleInternalForce(const ElementType & type, GhostType gt); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); /// TODO remove void computeRotationMatrix(const ElementType & type); protected: /// compute Rotation Matrices template void computeRotationMatrix(__attribute__((unused)) Array & rotations) {} /* ------------------------------------------------------------------------ */ /* Mass (structural_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// computes rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// finish the computation of residual to solve in increment void updateResidualInternal(); /* ------------------------------------------------------------------------ */ private: template void assembleStiffnessMatrix(); template void assembleMass(); template void computeStressOnQuad(); template void computeTangentModuli(Array & tangent_moduli); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: - dumper::Field * createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - dumper::Field * createNodalFieldBool(const std::string & field_name, - const std::string & group_name, - bool padding_flag) override; - - dumper::Field * createElementalField(const std::string & field_name, - const std::string & group_name, - bool padding_flag, - const UInt & spatial_dimension, - const ElementKind & kind) override; + std::shared_ptr + createNodalFieldReal(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; + + std::shared_ptr + createNodalFieldBool(const std::string & field_name, + const std::string & group_name, + bool padding_flag) override; + + std::shared_ptr + createElementalField(const std::string & field_name, + const std::string & group_name, bool padding_flag, + const UInt & spatial_dimension, + const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step // void setTimeStep(Real time_step, const ID & solver_id = "") override; /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the StructuralMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement_rotation, Array &); /// get the StructuralMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the StructuralMechanicsModel::acceleration vector, updated /// by /// StructuralMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the StructuralMechanicsModel::external_force vector AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the StructuralMechanicsModel::internal_force vector (boundary forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the StructuralMechanicsModel::boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(RotationMatrix, rotation_matrix, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Set_ID, set_ID, UInt); void addMaterial(StructuralMaterial & material) { materials.push_back(material); } const StructuralMaterial & getMaterial(const Element & element) const { return materials[element_material(element)]; } /* ------------------------------------------------------------------------ */ /* Boundaries (structural_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// Compute Linear load function set in global axis template void computeForcesByGlobalTractionArray(const Array & tractions); /// Compute Linear load function set in local axis template void computeForcesByLocalTractionArray(const Array & tractions); /// compute force vector from a function(x,y,momentum) that describe stresses // template // void computeForcesFromFunction(BoundaryFunction in_function, // BoundaryFunctionType function_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement_rotation{nullptr}; /// velocities array Array * velocity{nullptr}; /// accelerations array Array * acceleration{nullptr}; /// forces array Array * internal_force{nullptr}; /// forces array Array * external_force{nullptr}; /// lumped mass array Array * mass{nullptr}; /// boundaries array Array * blocked_dofs{nullptr}; /// stress array ElementTypeMapArray stress; ElementTypeMapArray element_material; // Define sets of beams ElementTypeMapArray set_ID; /// number of degre of freedom UInt nb_degree_of_freedom; // Rotation matrix ElementTypeMapArray rotation_matrix; // /// analysis method check the list in akantu::AnalysisMethod // AnalysisMethod method; /// flag defining if the increment must be computed or not bool increment_flag; /* ------------------------------------------------------------------------ */ std::vector materials; }; } // namespace akantu #endif /* __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ */ diff --git a/test/test_io/test_dumper/test_dumper.cc b/test/test_io/test_dumper/test_dumper.cc index 93a39389e..58e4d9c49 100644 --- a/test/test_io/test_dumper/test_dumper.cc +++ b/test/test_io/test_dumper/test_dumper.cc @@ -1,156 +1,157 @@ /** * @file test_dumper.cc * * @author David Simon Kammer * * @date creation: Tue Sep 02 2014 * @date last modification: Mon Jan 22 2018 * * @brief test dumper * * @section LICENSE * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" #include "dumper_variable.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize("input_file.dat", argc, argv); UInt spatial_dimension = 3; Mesh mesh(spatial_dimension); mesh.read("test_dumper.msh"); SolidMechanicsModel model(mesh); auto && mat_selector = std::make_shared>("physical_names", model); model.setMaterialSelector(mat_selector); model.initFull(); model.assembleInternalForces(); Real time_step = 0.1; const Array & coord = mesh.getNodes(); Array & disp = model.getDisplacement(); Array & bound = model.getBlockedDOFs(); for (UInt n = 0; n < mesh.getNbNodes(); ++n) { Real dist = 0.; for (UInt d = 0; d < spatial_dimension; ++d) { dist += coord(n, d) * coord(n, d); } dist = sqrt(dist); for (UInt d = 0; d < spatial_dimension; ++d) { disp(n, d) = (d + 1) * dist; bound(n, d) = bool((n % 2) + d); } } // dump boundary bottom as reference model.setGroupDirectory("paraview", "Bottom"); model.setGroupBaseName("paraview_bottom", "Bottom"); model.addDumpGroupField("displacement", "Bottom"); model.addDumpGroupField("blocked_dofs", "Bottom"); UInt nbp = 3; DumperParaview prvdumper("paraview_bottom_parallel", "paraview", false); iohelper::Dumper & prvdpr = prvdumper.getDumper(); for (UInt p = 0; p < nbp; ++p) { prvdpr.setParallelContext(p, nbp, 0); if (p != 0) { prvdumper.unRegisterField("connectivities"); prvdumper.unRegisterField("element_type"); prvdumper.unRegisterField("positions"); prvdumper.unRegisterField("displacement"); } prvdumper.registerFilteredMesh(mesh, mesh.getElementGroup("Bottom").getElements(), mesh.getElementGroup("Bottom").getNodes()); prvdumper.registerField("displacement", - new dumper::NodalField( + std::make_shared>( model.getDisplacement(), 0, 0, &(mesh.getElementGroup("Bottom").getNodes()))); prvdumper.dump(0); } DumperText txtdumper("text_bottom", iohelper::_tdm_csv); txtdumper.setDirectory("paraview"); txtdumper.setPrecision(8); txtdumper.setTimeStep(time_step); txtdumper.registerFilteredMesh(mesh, mesh.getElementGroup("Bottom").getElements(), mesh.getElementGroup("Bottom").getNodes()); txtdumper.registerField("displacement", - new dumper::NodalField( + std::make_shared>( model.getDisplacement(), 0, 0, &(mesh.getElementGroup("Bottom").getNodes()))); txtdumper.registerField("blocked_dofs", - new dumper::NodalField( + std::make_shared>( model.getBlockedDOFs(), 0, 0, &(mesh.getElementGroup("Bottom").getNodes()))); Real pot_energy = 1.2345567891; Vector gforces(2, 1.); - txtdumper.registerVariable("potential_energy", - new dumper::Variable(pot_energy)); - txtdumper.registerVariable("global_forces", - new dumper::Variable>(gforces)); + txtdumper.registerVariable( + "potential_energy", std::make_shared>(pot_energy)); + txtdumper.registerVariable( + "global_forces", + std::make_shared>>(gforces)); // dump a first time before the main loop model.dumpGroup("Bottom"); txtdumper.dump(); Real time = 0.; for (UInt i = 1; i < 5; ++i) { pot_energy += 2.; gforces(0) += 0.1; gforces(1) += 0.2; // pre -> cor // increment time after all steps of integration time += time_step; // dump after time increment if (i % 2 == 0) { txtdumper.dump(time, i); model.dumpGroup("Bottom"); // parallel test for (UInt p = 0; p < nbp; ++p) { prvdpr.setParallelContext(p, nbp, 0); prvdumper.dump(i); } } } finalize(); return EXIT_SUCCESS; } diff --git a/test/test_mesh_utils/test_mesh_partitionate/test_mesh_partitionate_scotch.cc b/test/test_mesh_utils/test_mesh_partitionate/test_mesh_partitionate_scotch.cc index 51acadcc3..506127980 100644 --- a/test/test_mesh_utils/test_mesh_partitionate/test_mesh_partitionate_scotch.cc +++ b/test/test_mesh_utils/test_mesh_partitionate/test_mesh_partitionate_scotch.cc @@ -1,75 +1,74 @@ /** * @file test_mesh_partitionate_scotch.cc * * @author Nicolas Richart * * @date creation: Sun Sep 12 2010 * @date last modification: Mon Jan 22 2018 * * @brief test of internal facet extraction * * @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_common.hh" #include "mesh.hh" #include "mesh_partition_scotch.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumper_elemental_field.hh" #include "dumper_iohelper_paraview.hh" #endif // AKANTU_USE_IOHELPER using namespace akantu; /* -------------------------------------------------------------------------- */ /* Main */ /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); debug::setDebugLevel(akantu::dblDump); int dim = 2; Mesh mesh(dim); mesh.read("triangle.msh"); MeshPartitionScotch partition(mesh, dim); partition.partitionate(8); #ifdef AKANTU_USE_IOHELPER DumperParaview dumper("test-scotch-partition"); - dumper::Field * field = - new dumper::ElementalField(partition.getPartitions(), dim); + auto field = std::make_shared>( + partition.getPartitions(), dim); dumper.registerMesh(mesh, dim); dumper.registerField("partitions", field); dumper.dump(); #endif // AKANTU_USE_IOHELPER partition.reorder(); mesh.write("triangle_reorder.msh"); - finalize(); return EXIT_SUCCESS; }