diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc index 6eebb30a7..f766d3624 100644 --- a/src/common/aka_error.cc +++ b/src/common/aka_error.cc @@ -1,355 +1,356 @@ /** * @file aka_error.cc * * @author Nicolas Richart * * @date creation: Mon Sep 06 2010 * @date last modification: Tue Jan 19 2016 * * @brief handling of errors * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_error.hh" #include "aka_common.hh" #include "aka_config.hh" /* -------------------------------------------------------------------------- */ #include #include #if (defined(READLINK_COMMAND) || defined(ADDR2LINE_COMMAND)) && \ (not defined(_WIN32)) #include #include #endif #include #include #include #include #include #include #include #include #if defined(AKANTU_CORE_CXX11) #include #elif defined(AKANTU_USE_OBSOLETE_GETTIMEOFDAY) #include #else #include #endif #ifdef AKANTU_USE_MPI #include #endif /* -------------------------------------------------------------------------- */ namespace akantu { namespace debug { static void printBacktraceAndExit(int) { std::terminate(); } /* ------------------------------------------------------------------------ */ void initSignalHandler() { std::signal(SIGSEGV, &printBacktraceAndExit); } /* ------------------------------------------------------------------------ */ std::string demangle(const char * symbol) { int status; std::string result; char * demangled_name; if ((demangled_name = abi::__cxa_demangle(symbol, nullptr, nullptr, &status)) != nullptr) { result = demangled_name; free(demangled_name); } else { result = symbol; } return result; } /* ------------------------------------------------------------------------ */ #if (defined(READLINK_COMMAND) || defined(ADDR2LINK_COMMAND)) && \ (not defined(_WIN32)) std::string exec(const std::string & cmd) { FILE * pipe = popen(cmd.c_str(), "r"); if (!pipe) return ""; char buffer[1024]; std::string result = ""; while (!feof(pipe)) { if (fgets(buffer, 128, pipe) != nullptr) result += buffer; } result = result.substr(0, result.size() - 1); pclose(pipe); return result; } #endif /* ------------------------------------------------------------------------ */ void printBacktrace(__attribute__((unused)) int sig) { AKANTU_DEBUG_INFO("Caught signal " << sig << "!"); #if not defined(_WIN32) #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) std::string me = ""; char buf[1024]; /* The manpage says it won't null terminate. Let's zero the buffer. */ memset(buf, 0, sizeof(buf)); /* Note we use sizeof(buf)-1 since we may need an extra char for NUL. */ if (readlink("/proc/self/exe", buf, sizeof(buf) - 1)) me = std::string(buf); std::ifstream inmaps; inmaps.open("/proc/self/maps"); std::map addr_map; std::string line; while (inmaps.good()) { std::getline(inmaps, line); std::stringstream sstr(line); size_t first = line.find('-'); std::stringstream sstra(line.substr(0, first)); size_t addr; sstra >> std::hex >> addr; std::string lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; if (lib != "" && addr_map.find(lib) == addr_map.end()) { addr_map[lib] = addr; } } if (me != "") addr_map[me] = 0; #endif /// \todo for windows this part could be coded using CaptureStackBackTrace /// and SymFromAddr const size_t max_depth = 100; size_t stack_depth; void * stack_addrs[max_depth]; char ** stack_strings; size_t i; stack_depth = backtrace(stack_addrs, max_depth); stack_strings = backtrace_symbols(stack_addrs, stack_depth); std::cerr << "BACKTRACE : " << stack_depth << " stack frames." << std::endl; auto w = size_t(std::floor(log(double(stack_depth)) / std::log(10.)) + 1); /// -1 to remove the call to the printBacktrace function for (i = 1; i < stack_depth; i++) { std::cerr << std::dec << " [" << std::setw(w) << i << "] "; std::string bt_line(stack_strings[i]); size_t first, second; if ((first = bt_line.find('(')) != std::string::npos && (second = bt_line.find('+')) != std::string::npos) { std::string location = bt_line.substr(0, first); #if defined(READLINK_COMMAND) - location = exec(std::string(BOOST_PP_STRINGIZE(READLINK_COMMAND)) + - std::string(" -f ") + location); + std::string location_cmd = std::string(BOOST_PP_STRINGIZE(READLINK_COMMAND)) + + std::string(" -f ") + location; + location = exec(location_cmd); #endif std::string call = demangle(bt_line.substr(first + 1, second - first - 1).c_str()); size_t f = bt_line.find('['); size_t s = bt_line.find(']'); std::string address = bt_line.substr(f + 1, s - f - 1); std::stringstream sstra(address); size_t addr; sstra >> std::hex >> addr; std::cerr << location << " [" << call << "]"; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) auto it = addr_map.find(location); if (it != addr_map.end()) { std::stringstream syscom; syscom << BOOST_PP_STRINGIZE(ADDR2LINE_COMMAND) << " 0x" << std::hex << (addr - it->second) << " -i -e " << location; std::string line = exec(syscom.str()); std::cerr << " (" << line << ")" << std::endl; } else { #endif std::cerr << " (0x" << std::hex << addr << ")" << std::endl; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) } #endif } else { std::cerr << bt_line << std::endl; } } free(stack_strings); std::cerr << "END BACKTRACE" << std::endl; #endif } /* ------------------------------------------------------------------------ */ namespace { void terminate_handler() { auto eptr = std::current_exception(); auto t = abi::__cxa_current_exception_type(); auto name = t ? demangle(t->name()) : std::string("unknown"); try { if (eptr) std::rethrow_exception(eptr); else std::cerr << "!! Execution terminated for unknown reasons !!" << std::endl; } catch (std::exception & e) { std::cerr << "!! Uncaught exception of type " << name << " !!\nwhat(): \"" << e.what() << "\"" << std::endl; } catch (...) { std::cerr << "!! Something strange of type \"" << name << "\" was thrown.... !!" << std::endl; } if (debugger.printBacktrace()) printBacktrace(15); } } // namespace /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ Debugger::Debugger() { cout = &std::cerr; level = dblWarning; parallel_context = ""; file_open = false; print_backtrace = false; initSignalHandler(); std::set_terminate(terminate_handler); } /* ------------------------------------------------------------------------ */ Debugger::~Debugger() { if (file_open) { dynamic_cast(cout)->close(); delete cout; } } /* ------------------------------------------------------------------------ */ void Debugger::exit(int status) { if (status != 0) std::terminate(); std::exit(0); } /*------------------------------------------------------------------------- */ void Debugger::throwException(const std::string & info, const std::string & file, unsigned int line, __attribute__((unused)) bool silent, __attribute__((unused)) const std::string & location) const noexcept(false) { #if !defined(AKANTU_NDEBUG) if (!silent) { printMessage("###", dblWarning, info + location); } #endif debug::Exception ex(info, file, line); throw ex; } /* ------------------------------------------------------------------------ */ void Debugger::printMessage(const std::string & prefix, const DebugLevel & level, const std::string & info) const { if (this->level >= level) { #if defined(AKANTU_CORE_CXX11) double timestamp = std::chrono::duration_cast>( std::chrono::system_clock::now().time_since_epoch()) .count(); #elif defined(AKANTU_USE_OBSOLETE_GETTIMEOFDAY) struct timeval time; gettimeofday(&time, NULL); double timestamp = time.tv_sec * 1e6 + time.tv_usec; /*in us*/ #else struct timespec time; clock_gettime(CLOCK_REALTIME_COARSE, &time); double timestamp = time.tv_sec * 1e6 + time.tv_nsec * 1e-3; /*in us*/ #endif *(cout) << parallel_context << "{" << (size_t)timestamp << "} " << prefix << " " << info << std::endl; } } /* ------------------------------------------------------------------------ */ void Debugger::setDebugLevel(const DebugLevel & level) { this->level = level; } /* ------------------------------------------------------------------------ */ const DebugLevel & Debugger::getDebugLevel() const { return this->level; } /* ------------------------------------------------------------------------ */ void Debugger::setLogFile(const std::string & filename) { if (file_open) { dynamic_cast(cout)->close(); delete cout; } auto * fileout = new std::ofstream(filename.c_str()); file_open = true; cout = fileout; } std::ostream & Debugger::getOutputStream() { return *cout; } /* ------------------------------------------------------------------------ */ void Debugger::setParallelContext(int rank, int size) { std::stringstream sstr; UInt pad = std::ceil(std::log10(size)); sstr << "<" << getpid() << ">[R" << std::setfill(' ') << std::right << std::setw(pad) << rank << "|S" << size << "] "; parallel_context = sstr.str(); } void setDebugLevel(const DebugLevel & level) { debugger.setDebugLevel(level); } const DebugLevel & getDebugLevel() { return debugger.getDebugLevel(); } /* -------------------------------------------------------------------------- */ void exit(int status) { debugger.exit(status); } } // namespace debug } // namespace akantu diff --git a/src/io/dumper/dumper_element_partition.hh b/src/io/dumper/dumper_element_partition.hh index f037a7c4a..d735814cb 100644 --- a/src/io/dumper/dumper_element_partition.hh +++ b/src/io/dumper/dumper_element_partition.hh @@ -1,123 +1,123 @@ /** * @file dumper_element_partition.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Sep 02 2014 * @date last modification: Tue Jul 14 2015 * * @brief ElementPartition field * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ namespace akantu { __BEGIN_AKANTU_DUMPER__ #ifdef AKANTU_IGFEM # include "dumper_igfem_element_partition.hh" #endif /* -------------------------------------------------------------------------- */ template class element_partition_field_iterator : public element_iterator { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ public: - typedef element_iterator parent; + using parent = element_iterator; using return_type = typename SingleType::return_type; using array_iterator = typename types::array_iterator; using field_type = typename types::field_type; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: element_partition_field_iterator(const field_type & field, const typename field_type::type_iterator & t_it, const typename field_type::type_iterator & t_it_end, const array_iterator & array_it, const array_iterator & array_it_end, const GhostType ghost_type = _not_ghost) : parent(field, t_it, t_it_end, array_it, array_it_end, ghost_type) { prank = Communicator::getStaticCommunicator().whoAmI(); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: return_type operator*() { return return_type(1, prank); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: UInt prank; }; /* -------------------------------------------------------------------------- */ template class ElementPartitionField : public GenericElementalField, element_partition_field_iterator> { public: /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ - typedef SingleType types; + using types = SingleType; using iterator = element_partition_field_iterator; - typedef GenericElementalField parent; + using parent = GenericElementalField; using field_type = typename types::field_type; public: /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ ElementPartitionField(const field_type & field, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) : parent(field, spatial_dimension, ghost_type, element_kind) { this->homogeneous = true; } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ UInt getDim() override { return 1; } }; /* -------------------------------------------------------------------------- */ __END_AKANTU_DUMPER__ } // akantu diff --git a/src/io/mesh_io/mesh_io_diana.cc b/src/io/mesh_io/mesh_io_diana.cc index 0e2ec29b6..089edfdde 100644 --- a/src/io/mesh_io/mesh_io_diana.cc +++ b/src/io/mesh_io/mesh_io_diana.cc @@ -1,602 +1,602 @@ /** * @file mesh_io_diana.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Nicolas Richart * @author Alodie Schneuwly * * @date creation: Sat Mar 26 2011 * @date last modification: Thu Jan 21 2016 * * @brief handles diana meshes * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "element_group.hh" #include "mesh_io_diana.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include namespace akantu { /* -------------------------------------------------------------------------- */ /* Methods Implentations */ /* -------------------------------------------------------------------------- */ MeshIODiana::MeshIODiana() { canReadSurface = true; canReadExtendedData = true; _diana_to_akantu_element_types["T9TM"] = _triangle_3; _diana_to_akantu_element_types["CT6CM"] = _triangle_6; _diana_to_akantu_element_types["Q12TM"] = _quadrangle_4; _diana_to_akantu_element_types["CQ8CM"] = _quadrangle_8; _diana_to_akantu_element_types["TP18L"] = _pentahedron_6; _diana_to_akantu_element_types["CTP45"] = _pentahedron_15; _diana_to_akantu_element_types["TE12L"] = _tetrahedron_4; _diana_to_akantu_element_types["HX24L"] = _hexahedron_8; _diana_to_akantu_element_types["CHX60"] = _hexahedron_20; _diana_to_akantu_mat_prop["YOUNG"] = "E"; _diana_to_akantu_mat_prop["DENSIT"] = "rho"; _diana_to_akantu_mat_prop["POISON"] = "nu"; std::map::iterator it; for (it = _diana_to_akantu_element_types.begin(); it != _diana_to_akantu_element_types.end(); ++it) { UInt nb_nodes = Mesh::getNbNodesPerElement(it->second); auto * tmp = new UInt[nb_nodes]; for (UInt i = 0; i < nb_nodes; ++i) { tmp[i] = i; } switch (it->second) { case _tetrahedron_10: tmp[8] = 9; tmp[9] = 8; break; case _pentahedron_15: tmp[0] = 2; tmp[1] = 8; tmp[2] = 0; tmp[3] = 6; tmp[4] = 1; tmp[5] = 7; tmp[6] = 11; tmp[7] = 9; tmp[8] = 10; tmp[9] = 5; tmp[10] = 14; tmp[11] = 3; tmp[12] = 12; tmp[13] = 4; tmp[14] = 13; break; case _hexahedron_20: tmp[0] = 5; tmp[1] = 16; tmp[2] = 4; tmp[3] = 19; tmp[4] = 7; tmp[5] = 18; tmp[6] = 6; tmp[7] = 17; tmp[8] = 13; tmp[9] = 12; tmp[10] = 15; tmp[11] = 14; tmp[12] = 1; tmp[13] = 8; tmp[14] = 0; tmp[15] = 11; tmp[16] = 3; tmp[17] = 10; tmp[18] = 2; tmp[19] = 9; break; default: // nothing to change break; } _read_order[it->second] = tmp; } } /* -------------------------------------------------------------------------- */ MeshIODiana::~MeshIODiana() = default; /* -------------------------------------------------------------------------- */ inline void my_getline(std::ifstream & infile, std::string & line) { std::getline(infile, line); // read the line size_t pos = line.find('\r'); /// remove the extra \r if needed line = line.substr(0, pos); } /* -------------------------------------------------------------------------- */ void MeshIODiana::read(const std::string & filename, Mesh & mesh) { AKANTU_DEBUG_IN(); MeshAccessor mesh_accessor(mesh); std::ifstream infile; infile.open(filename.c_str()); std::string line; UInt first_node_number = std::numeric_limits::max(); diana_element_number_to_elements.clear(); if (!infile.good()) { AKANTU_DEBUG_ERROR("Cannot open file " << filename); } while (infile.good()) { my_getline(infile, line); /// read all nodes if (line == "'COORDINATES'") { line = readCoordinates(infile, mesh, first_node_number); } /// read all elements if (line == "'ELEMENTS'") { line = readElements(infile, mesh, first_node_number); } /// read the material properties and write a .dat file if (line == "'MATERIALS'") { line = readMaterial(infile, filename); } /// read the material properties and write a .dat file if (line == "'GROUPS'") { line = readGroups(infile, mesh, first_node_number); } } infile.close(); mesh_accessor.setNbGlobalNodes(mesh.getNbNodes()); MeshUtils::fillElementToSubElementsData(mesh); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshIODiana::write(__attribute__((unused)) const std::string & filename, __attribute__((unused)) const Mesh & mesh) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readCoordinates(std::ifstream & infile, Mesh & mesh, UInt & first_node_number) { AKANTU_DEBUG_IN(); MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); std::string line; UInt index; Vector coord(3); do { my_getline(infile, line); if ("'ELEMENTS'" == line) break; std::stringstream sstr_node(line); sstr_node >> index >> coord(0) >> coord(1) >> coord(2); first_node_number = first_node_number < index ? first_node_number : index; nodes.push_back(coord); } while (true); AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ UInt MeshIODiana::readInterval(std::stringstream & line, std::set & interval) { UInt first; line >> first; if (line.fail()) { return 0; } interval.insert(first); UInt second; int dash; dash = line.get(); if (dash == '-') { line >> second; interval.insert(second); return 2; } if (line.fail()) line.clear(std::ios::eofbit); // in case of get at end of the line else line.unget(); return 1; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readGroups(std::ifstream & infile, Mesh & mesh, UInt first_node_number) { AKANTU_DEBUG_IN(); std::string line; my_getline(infile, line); bool reading_nodes_group = false; while (line != "'SUPPORTS'") { if (line == "NODES") { reading_nodes_group = true; my_getline(infile, line); } if (line == "ELEMEN") { reading_nodes_group = false; my_getline(infile, line); } auto * str = new std::stringstream(line); UInt id; std::string name; char c; *str >> id >> name >> c; auto * list_ids = new Array(0, 1, name); UInt s = 1; bool end = false; while (!end) { while (!str->eof() && s != 0) { std::set interval; s = readInterval(*str, interval); auto it = interval.begin(); if (s == 1) list_ids->push_back(*it); if (s == 2) { UInt first = *it; ++it; UInt second = *it; for (UInt i = first; i <= second; ++i) { list_ids->push_back(i); } } } if (str->fail()) end = true; else { my_getline(infile, line); delete str; str = new std::stringstream(line); } } delete str; if (reading_nodes_group) { NodeGroup & ng = mesh.createNodeGroup(name); for (UInt i = 0; i < list_ids->size(); ++i) { UInt node = (*list_ids)(i)-first_node_number; ng.add(node, false); } delete list_ids; } else { ElementGroup & eg = mesh.createElementGroup(name); for (UInt i = 0; i < list_ids->size(); ++i) { Element & elem = diana_element_number_to_elements[(*list_ids)(i)]; if (elem.type != _not_defined) eg.add(elem, false, false); } eg.optimize(); delete list_ids; } my_getline(infile, line); } AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readElements(std::ifstream & infile, Mesh & mesh, UInt first_node_number) { AKANTU_DEBUG_IN(); std::string line; my_getline(infile, line); if ("CONNECTIVITY" == line) { line = readConnectivity(infile, mesh, first_node_number); } /// read the line corresponding to the materials if ("MATERIALS" == line) { line = readMaterialElement(infile, mesh); } AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readConnectivity(std::ifstream & infile, Mesh & mesh, UInt first_node_number) { AKANTU_DEBUG_IN(); MeshAccessor mesh_accessor(mesh); Int index; std::string lline; std::string diana_type; ElementType akantu_type, akantu_type_old = _not_defined; Array * connectivity = nullptr; UInt node_per_element = 0; Element elem; UInt * read_order = nullptr; while (true) { my_getline(infile, lline); // std::cerr << lline << std::endl; std::stringstream sstr_elem(lline); if (lline == "MATERIALS") break; /// traiter les coordonnees sstr_elem >> index; sstr_elem >> diana_type; akantu_type = _diana_to_akantu_element_types[diana_type]; if (akantu_type == _not_defined) continue; if (akantu_type != akantu_type_old) { connectivity = &(mesh_accessor.getConnectivity(akantu_type)); node_per_element = connectivity->getNbComponent(); akantu_type_old = akantu_type; read_order = _read_order[akantu_type]; } Vector local_connect(node_per_element); // used if element is written on two lines UInt j_last = 0; for (UInt j = 0; j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; // check s'il y a pas plus rien après un certain point if (sstr_elem.fail()) { sstr_elem.clear(); sstr_elem.ignore(); break; } node_index -= first_node_number; local_connect(read_order[j]) = node_index; j_last = j; } // check if element is written in two lines if (j_last != (node_per_element - 1)) { // if this is the case, read on more line my_getline(infile, lline); std::stringstream sstr_elem(lline); for (UInt j = (j_last + 1); j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; node_index -= first_node_number; local_connect(read_order[j]) = node_index; } } connectivity->push_back(local_connect); elem.type = akantu_type; elem.element = connectivity->size() - 1; diana_element_number_to_elements[index] = elem; akantu_number_to_diana_number[elem] = index; } AKANTU_DEBUG_OUT(); return lline; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readMaterialElement(std::ifstream & infile, Mesh & mesh) { AKANTU_DEBUG_IN(); std::string line; std::stringstream sstr_tag_name; sstr_tag_name << "tag_" << 0; Mesh::type_iterator it = mesh.firstType(); Mesh::type_iterator end = mesh.lastType(); for (; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); mesh.getDataPointer("material", *it, _not_ghost, 1) .resize(nb_element); } my_getline(infile, line); while (line != "'MATERIALS'") { line = line.substr(line.find('/') + 1, std::string::npos); // erase the first slash / of the line char tutu[250]; strncpy(tutu, line.c_str(), 250); AKANTU_DEBUG_WARNING("reading line " << line); Array temp_id(0, 2); UInt mat; while (true) { std::stringstream sstr_intervals_elements(line); Vector id(2); char temp; while (sstr_intervals_elements.good()) { sstr_intervals_elements >> id(0) >> temp >> id(1); // >> "/" >> mat; if (!sstr_intervals_elements.fail()) temp_id.push_back(id); } if (sstr_intervals_elements.fail()) { sstr_intervals_elements.clear(); sstr_intervals_elements.ignore(); sstr_intervals_elements >> mat; break; } my_getline(infile, line); } // loop over elements // UInt * temp_id_val = temp_id.storage(); for (UInt i = 0; i < temp_id.size(); ++i) for (UInt j = temp_id(i, 0); j <= temp_id(i, 1); ++j) { Element & element = diana_element_number_to_elements[j]; if (element.type == _not_defined) continue; UInt elem = element.element; ElementType type = element.type; Array & data = mesh.getDataPointer("material", type, _not_ghost); data(elem) = mat; } my_getline(infile, line); } AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ std::string MeshIODiana::readMaterial(std::ifstream & infile, const std::string & filename) { AKANTU_DEBUG_IN(); std::stringstream mat_file_name; mat_file_name << "material_" << filename; std::ofstream material_file; material_file.open(mat_file_name.str().c_str()); // mat_file_name.str()); UInt mat_index; std::string line; bool first_mat = true; bool end = false; UInt mat_id = 0; - typedef std::map MatProp; + using MatProp = std::map; MatProp mat_prop; do { my_getline(infile, line); std::stringstream sstr_material(line); if (("'GROUPS'" == line) || ("'END'" == line)) { if (!mat_prop.empty()) { material_file << "material elastic [" << std::endl; material_file << "\tname = material" << ++mat_id << std::endl; for (auto it = mat_prop.begin(); it != mat_prop.end(); ++it) material_file << "\t" << it->first << " = " << it->second << std::endl; material_file << "]" << std::endl; mat_prop.clear(); } end = true; } else { /// traiter les caractéristiques des matériaux sstr_material >> mat_index; if (!sstr_material.fail()) { if (!first_mat) { if (!mat_prop.empty()) { material_file << "material elastic [" << std::endl; material_file << "\tname = material" << ++mat_id << std::endl; for (auto it = mat_prop.begin(); it != mat_prop.end(); ++it) material_file << "\t" << it->first << " = " << it->second << std::endl; material_file << "]" << std::endl; mat_prop.clear(); } } first_mat = false; } else { sstr_material.clear(); } std::string prop_name; sstr_material >> prop_name; std::map::iterator it; it = _diana_to_akantu_mat_prop.find(prop_name); if (it != _diana_to_akantu_mat_prop.end()) { Real value; sstr_material >> value; mat_prop[it->second] = value; } else { AKANTU_DEBUG_INFO("In material reader, property " << prop_name << "not recognized"); } } } while (!end); AKANTU_DEBUG_OUT(); return line; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/io/parser/parser_grammar_tmpl.hh b/src/io/parser/parser_grammar_tmpl.hh index fc04f6418..fa7fd7191 100644 --- a/src/io/parser/parser_grammar_tmpl.hh +++ b/src/io/parser/parser_grammar_tmpl.hh @@ -1,81 +1,81 @@ /** * @file parser_grammar_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Nov 11 2015 * * @brief implementation of the templated part of ParsableParam Parsable and * ParsableParamTyped * * @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 #include #include #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_PARSER_GRAMMAR_TMPL_HH #define AKANTU_PARSER_GRAMMAR_TMPL_HH namespace akantu { namespace qi = boost::spirit::qi; /* -------------------------------------------------------------------------- */ template T Parser::parseType(const std::string & value, Grammar & grammar) { using boost::spirit::ascii::space; std::string::const_iterator b = value.begin(); std::string::const_iterator e = value.end(); T resultat = T(); bool res = false; try { res = qi::phrase_parse(b, e, grammar, space, resultat); } catch(debug::Exception & ex) { AKANTU_EXCEPTION("Could not parse '" << value << "' as a " << debug::demangle(typeid(T).name()) << ", an unknown error append '" << ex.what()); } if(!res || (b != e)) { AKANTU_EXCEPTION("Could not parse '" << value << "' as a " << debug::demangle(typeid(T).name()) << ", an unknown error append '" << std::string(value.begin(), b) << "" << std::string(b, e) << "'"); } return resultat; } namespace parser { template struct Skipper { - typedef qi::rule type; + using type = qi::rule; }; } } // akantu #endif //AKANTU_PARSER_GRAMMAR_TMPL_HH diff --git a/src/mesh/group_manager_inline_impl.cc b/src/mesh/group_manager_inline_impl.cc index ad0170672..437abda2b 100644 --- a/src/mesh/group_manager_inline_impl.cc +++ b/src/mesh/group_manager_inline_impl.cc @@ -1,207 +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: Tue Dec 08 2015 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dumper_field.hh" #include "element_group.hh" #include "element_type_map_filter.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_nodal_field.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template class dump_type> dumper::Field * GroupManager::createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const ElementTypeMapArray * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") - return this->createElementalField >( + return this->createElementalField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); else - return this->createElementalFilteredField >( + return this->createElementalFilteredField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); } /* -------------------------------------------------------------------------- */ template class T2, template class, bool> class dump_type> dumper::Field * GroupManager::createElementalField( const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const ElementTypeMapArray * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") - return this->createElementalField >( + return this->createElementalField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); else - return this->createElementalFilteredField >( + return this->createElementalFilteredField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); } /* -------------------------------------------------------------------------- */ template class dump_type> ///< type of InternalMaterialField dumper::Field * GroupManager::createElementalField(const ElementTypeMapArray & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const ElementTypeMapArray * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") - return this->createElementalField >( + return this->createElementalField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); else - return this->createElementalFilteredField >( + return this->createElementalFilteredField>( field, group_name, spatial_dimension, kind, nb_data_per_elem); } /* -------------------------------------------------------------------------- */ template dumper::Field * GroupManager::createElementalField( const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, const ElementTypeMap & nb_data_per_elem) { const field_type * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name != "all") throw; dumper::Field * dumper = new dump_type(field, spatial_dimension, _not_ghost, kind); dumper->setNbDataPerElem(nb_data_per_elem); return dumper; } /* -------------------------------------------------------------------------- */ template dumper::Field * GroupManager::createElementalFilteredField( const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap nb_data_per_elem) { const field_type * field_ptr = &field; if (field_ptr == nullptr) return nullptr; if (group_name == "all") throw; using T = typename field_type::type; ElementGroup & group = this->getElementGroup(group_name); UInt dim = group.getDimension(); if (dim != spatial_dimension) throw; const ElementTypeMapArray & elemental_filter = group.getElements(); auto * filtered = new ElementTypeMapArrayFilter(field, elemental_filter, nb_data_per_elem); dumper::Field * dumper = new dump_type(*filtered, dim, _not_ghost, kind); dumper->setNbDataPerElem(nb_data_per_elem); return dumper; } /* -------------------------------------------------------------------------- */ - template class ftype> dumper::Field * GroupManager::createNodalField(const ftype * field, const std::string & group_name, UInt padding_size) { if (field == nullptr) return nullptr; if (group_name == "all") { - typedef typename dumper::NodalField DumpType; + using DumpType = typename dumper::NodalField; auto * dumper = new DumpType(*field, 0, 0, nullptr); dumper->setPadding(padding_size); return dumper; } else { ElementGroup & group = this->getElementGroup(group_name); const Array * nodal_filter = &(group.getNodes()); - typedef typename dumper::NodalField DumpType; + using DumpType = typename dumper::NodalField; auto * dumper = new DumpType(*field, 0, 0, nodal_filter); dumper->setPadding(padding_size); return dumper; } return nullptr; } /* -------------------------------------------------------------------------- */ - template class ftype> dumper::Field * GroupManager::createStridedNodalField(const ftype * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size) { if (field == NULL) return nullptr; if (group_name == "all") { - typedef typename dumper::NodalField DumpType; + using DumpType= typename dumper::NodalField; auto * dumper = new DumpType(*field, size, stride, NULL); dumper->setPadding(padding_size); return dumper; } else { ElementGroup & group = this->getElementGroup(group_name); const Array * nodal_filter = &(group.getNodes()); - typedef typename dumper::NodalField DumpType; + using DumpType = typename dumper::NodalField ; auto * dumper = new DumpType(*field, size, stride, nodal_filter); dumper->setPadding(padding_size); return dumper; } return nullptr; } /* -------------------------------------------------------------------------- */ -} // akantu +} // namespace akantu #endif diff --git a/src/mesh_utils/mesh_partition.cc b/src/mesh_utils/mesh_partition.cc index 4f73ff43e..b027c97e9 100644 --- a/src/mesh_utils/mesh_partition.cc +++ b/src/mesh_utils/mesh_partition.cc @@ -1,415 +1,415 @@ /** * @file mesh_partition.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Tue Aug 17 2010 * @date last modification: Fri Jan 22 2016 * * @brief implementation of common part of all partitioner * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_partition.hh" #include "aka_iterators.hh" #include "aka_types.hh" #include "mesh_accessor.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MeshPartition::MeshPartition(const Mesh & mesh, UInt spatial_dimension, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), mesh(mesh), spatial_dimension(spatial_dimension), partitions("partition", id, memory_id), ghost_partitions("ghost_partition", id, memory_id), ghost_partitions_offset("ghost_partition_offset", id, memory_id), saved_connectivity("saved_connectivity", id, memory_id), lin_to_element(mesh.getNbElement(spatial_dimension)) { AKANTU_DEBUG_IN(); Element elem; elem.ghost_type = _not_ghost; lin_to_element.resize(0); for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { elem.type = type; for (auto && i : arange(mesh.getNbElement(elem.type))) { elem.element = i; lin_to_element.push_back(elem); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MeshPartition::~MeshPartition() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * TODO this function should most probably be rewritten in a more modern way * conversion in c++ of the GENDUALMETIS (mesh.c) function wrote by George in * Metis (University of Minnesota) */ void MeshPartition::buildDualGraph(Array & dxadj, Array & dadjncy, Array & edge_loads, const EdgeLoadFunctor & edge_load_func) { AKANTU_DEBUG_IN(); // tweak mesh; UInt nb_good_types = 0; std::vector nb_nodes_per_element_p1; std::vector magic_number; std::vector nb_element; std::vector *> connectivities; UInt spatial_dimension = mesh.getSpatialDimension(); for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { ElementType type_p1 = Mesh::getP1ElementType(type); connectivities.push_back( &const_cast &>(mesh.getConnectivity(type, _not_ghost))); nb_nodes_per_element_p1.push_back(Mesh::getNbNodesPerElement(type_p1)); nb_element.push_back(connectivities.back()->size()); magic_number.push_back( Mesh::getNbNodesPerElement(Mesh::getFacetType(type_p1))); ++nb_good_types; } CSR node_to_elem; MeshUtils::buildNode2Elements(mesh, node_to_elem); UInt nb_total_element = std::accumulate(nb_element.begin(), nb_element.end(), 0); dxadj.resize(nb_total_element + 1); /// initialize the dxadj array auto dxadj_it = dxadj.begin(); for (auto && t : arange(nb_good_types)) { std::fill_n(dxadj_it, nb_element[t], nb_nodes_per_element_p1[t]); dxadj_it += nb_element[t]; } /// convert the dxadj_val array in a csr one for (UInt i = 1; i < nb_total_element; ++i) dxadj(i) += dxadj(i - 1); for (UInt i = nb_total_element; i > 0; --i) dxadj(i) = dxadj(i - 1); dxadj(0) = 0; dadjncy.resize(2 * dxadj(nb_total_element)); edge_loads.resize(2 * dxadj(nb_total_element)); /// weight map to determine adjacency std::unordered_map weight_map; for (UInt t = 0, linerized_el = 0; t < nb_good_types; ++t) { auto conn_it = connectivities[t]->begin(connectivities[t]->getNbComponent()); auto conn_end = connectivities[t]->end(connectivities[t]->getNbComponent()); for (; conn_it != conn_end; ++conn_it, ++linerized_el) { /// fill the weight map const auto & conn = *conn_it; for (UInt n : arange(nb_nodes_per_element_p1[t])) { auto && node = conn(n); for (auto k = node_to_elem.rbegin(node); k != node_to_elem.rend(node); --k) { auto && current_element = *k; auto && current_el = lin_to_element.find(current_element); if (current_el <= linerized_el) break; auto && weight_map_insert = weight_map.insert(std::make_pair(current_el, 1)); if (not weight_map_insert.second) (weight_map_insert.first->second)++; } } /// each element with a weight of the size of a facet are adjacent for (auto && weight_pair : weight_map) { UInt magic, adjacent_el; std::tie(adjacent_el, magic) = weight_pair; if (magic == magic_number[t]) { #if defined(AKANTU_COHESIVE_ELEMENT) /// Patch in order to prevent neighboring cohesive elements /// from detecting each other ElementKind linearized_el_kind = lin_to_element(linerized_el).kind(); ElementKind adjacent_el_kind = lin_to_element(adjacent_el).kind(); if (linearized_el_kind == adjacent_el_kind && linearized_el_kind == _ek_cohesive) continue; #endif UInt index_adj = dxadj(adjacent_el)++; UInt index_lin = dxadj(linerized_el)++; dadjncy(index_lin) = adjacent_el; dadjncy(index_adj) = linerized_el; } } weight_map.clear(); } } Int k_start = 0; for (UInt t = 0, linerized_el = 0, j = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) { for (Int k = k_start; k < dxadj(linerized_el); ++k, ++j) dadjncy(j) = dadjncy(k); dxadj(linerized_el) = j; k_start += nb_nodes_per_element_p1[t]; } for (UInt i = nb_total_element; i > 0; --i) dxadj(i) = dxadj(i - 1); dxadj(0) = 0; UInt adj = 0; for (UInt i = 0; i < nb_total_element; ++i) { UInt nb_adj = dxadj(i + 1) - dxadj(i); for (UInt j = 0; j < nb_adj; ++j, ++adj) { Int el_adj_id = dadjncy(dxadj(i) + j); Element el = lin_to_element(i); Element el_adj = lin_to_element(el_adj_id); Int load = edge_load_func(el, el_adj); edge_loads(adj) = load; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartition::fillPartitionInformation( const Mesh & mesh, const Int * linearized_partitions) { AKANTU_DEBUG_IN(); CSR node_to_elem; MeshUtils::buildNode2Elements(mesh, node_to_elem); UInt linearized_el = 0; for (auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); partitions.alloc(nb_element, 1, type, _not_ghost); auto & ghost_part_csr = ghost_partitions_csr(type, _not_ghost); ghost_part_csr.resizeRows(nb_element); ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost); ghost_partitions.alloc(0, 1, type, _ghost); const Array & connectivity = mesh.getConnectivity(type, _not_ghost); for (UInt el = 0; el < nb_element; ++el, ++linearized_el) { UInt part = linearized_partitions[linearized_el]; partitions(type, _not_ghost)(el) = part; std::list list_adj_part; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity.storage()[el * nb_nodes_per_element + n]; CSR::iterator ne; for (ne = node_to_elem.begin(node); ne != node_to_elem.end(node); ++ne) { const Element & adj_element = *ne; UInt adj_el = lin_to_element.find(adj_element); UInt adj_part = linearized_partitions[adj_el]; if (part != adj_part) { list_adj_part.push_back(adj_part); } } } list_adj_part.sort(); list_adj_part.unique(); for (auto & adj_part : list_adj_part) { ghost_part_csr.getRows().push_back(adj_part); ghost_part_csr.rowOffset(el)++; ghost_partitions(type, _ghost).push_back(adj_part); ghost_partitions_offset(type, _ghost)(el)++; } } ghost_part_csr.countToCSR(); /// convert the ghost_partitions_offset array in an offset array Array & ghost_partitions_offset_ptr = ghost_partitions_offset(type, _ghost); for (UInt i = 1; i < nb_element; ++i) ghost_partitions_offset_ptr(i) += ghost_partitions_offset_ptr(i - 1); for (UInt i = nb_element; i > 0; --i) ghost_partitions_offset_ptr(i) = ghost_partitions_offset_ptr(i - 1); ghost_partitions_offset_ptr(0) = 0; } // All Facets for (Int sp = spatial_dimension - 1; sp >= 0; --sp) { for (auto & type : mesh.elementTypes(sp, _not_ghost, _ek_not_defined)) { UInt nb_element = mesh.getNbElement(type); partitions.alloc(nb_element, 1, type, _not_ghost); AKANTU_DEBUG_INFO("Allocating partitions for " << type); auto & ghost_part_csr = ghost_partitions_csr(type, _not_ghost); ghost_part_csr.resizeRows(nb_element); ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost); ghost_partitions.alloc(0, 1, type, _ghost); AKANTU_DEBUG_INFO("Allocating ghost_partitions for " << type); const Array> & elem_to_subelem = mesh.getElementToSubelement(type, _not_ghost); for (UInt i(0); i < mesh.getNbElement(type, _not_ghost); ++i) { // Facet loop const std::vector & adjacent_elems = elem_to_subelem(i); if (!adjacent_elems.empty()) { - Element min_elem; + Element min_elem{_max_element_type, std::numeric_limits::max(), *ghost_type_t::end()}; UInt min_part(std::numeric_limits::max()); std::set adjacent_parts; for (UInt j(0); j < adjacent_elems.size(); ++j) { UInt adjacent_elem_id = adjacent_elems[j].element; UInt adjacent_elem_part = partitions(adjacent_elems[j].type, adjacent_elems[j].ghost_type)(adjacent_elem_id); if (adjacent_elem_part < min_part) { min_part = adjacent_elem_part; min_elem = adjacent_elems[j]; } adjacent_parts.insert(adjacent_elem_part); } partitions(type, _not_ghost)(i) = min_part; auto git = ghost_partitions_csr(min_elem.type, _not_ghost) .begin(min_elem.element); auto gend = ghost_partitions_csr(min_elem.type, _not_ghost) .end(min_elem.element); for (; git != gend; ++git) { adjacent_parts.insert(min_part); } adjacent_parts.erase(min_part); for (auto & part : adjacent_parts) { ghost_part_csr.getRows().push_back(part); ghost_part_csr.rowOffset(i)++; ghost_partitions(type, _ghost).push_back(part); } ghost_partitions_offset(type, _ghost)(i + 1) = ghost_partitions_offset(type, _ghost)(i + 1) + adjacent_elems.size(); } else { partitions(type, _not_ghost)(i) = 0; } } ghost_part_csr.countToCSR(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartition::tweakConnectivity(const Array & pairs) { AKANTU_DEBUG_IN(); if (pairs.size() == 0) return; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; Array & conn = const_cast &>(mesh.getConnectivity(type, _not_ghost)); UInt nb_nodes_per_element = conn.getNbComponent(); UInt nb_element = conn.size(); Array & saved_conn = saved_connectivity.alloc( nb_element, nb_nodes_per_element, type, _not_ghost); saved_conn.copy(conn); for (UInt i = 0; i < pairs.size(); ++i) { for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { if (pairs(i, 1) == conn(el, n)) conn(el, n) = pairs(i, 0); } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartition::restoreConnectivity() { AKANTU_DEBUG_IN(); MeshAccessor mesh_accessor(const_cast(mesh)); for (auto && type : saved_connectivity.elementTypes( spatial_dimension, _not_ghost, _ek_not_defined)) { auto & conn = mesh_accessor.getConnectivity(type, _not_ghost); auto & saved_conn = saved_connectivity(type, _not_ghost); conn.copy(saved_conn); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool MeshPartition::hasPartitions(const ElementType & type, const GhostType & ghost_type) { return partitions.exists(type, ghost_type); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index c08614023..7f29473d2 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1588 +1,1588 @@ /** * @file material.cc * * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Nov 24 2015 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(_st_material, id), is_init(false), fem(model.getFEEngine()), finite_deformation(false), name(""), model(model), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id, this->memory_id), stress("stress", *this), eigengradu("eigen_grad_u", *this), gradu("grad_u", *this), green_strain("green_strain", *this), piola_kirchhoff_2("piola_kirchhoff_2", *this), potential_energy("potential_energy", *this), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse coordinates", *this), interpolation_points_matrices("interpolation points matrices", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(model.getMesh(), _spatial_dimension = spatial_dimension); // model.getMesh().initElementTypeMapArray(element_filter, 1, // spatial_dimension, // false, _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Memory(id, model.getMemoryID()), Parsable(_st_material, id), is_init(false), fem(model.getFEEngine()), finite_deformation(false), name(""), model(model), spatial_dimension(dim), element_filter("element_filter", id, this->memory_id), stress("stress", *this, dim, fe_engine, this->element_filter), eigengradu("eigen_grad_u", *this, dim, fe_engine, this->element_filter), gradu("gradu", *this, dim, fe_engine, this->element_filter), green_strain("green_strain", *this, dim, fe_engine, this->element_filter), piola_kirchhoff_2("poila_kirchhoff_2", *this, dim, fe_engine, this->element_filter), potential_energy("potential_energy", *this, dim, fe_engine, this->element_filter), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse_coordinates", *this, dim, fe_engine, this->element_filter), interpolation_points_matrices("interpolation points matrices", *this, dim, fe_engine, this->element_filter) { AKANTU_DEBUG_IN(); element_filter.initialize(mesh, _spatial_dimension = spatial_dimension); // mesh.initElementTypeMapArray(element_filter, 1, spatial_dimension, false, // _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() = default; /* -------------------------------------------------------------------------- */ void Material::initialize() { registerParam("rho", rho, Real(0.), _pat_parsable | _pat_modifiable, "Density"); registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("finite_deformation", finite_deformation, false, _pat_parsable | _pat_readable, "Is finite deformation"); registerParam("inelastic_deformation", inelastic_deformation, false, _pat_internal, "Is inelastic deformation"); /// allocate gradu stress for local elements eigengradu.initialize(spatial_dimension * spatial_dimension); gradu.initialize(spatial_dimension * spatial_dimension); stress.initialize(spatial_dimension * spatial_dimension); potential_energy.initialize(1); this->model.registerEventHandler(*this); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); if (finite_deformation) { this->piola_kirchhoff_2.initialize(spatial_dimension * spatial_dimension); if (use_previous_stress) this->piola_kirchhoff_2.initializeHistory(); this->green_strain.initialize(spatial_dimension * spatial_dimension); } if (use_previous_stress) this->stress.initializeHistory(); if (use_previous_gradu) this->gradu.initializeHistory(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::savePreviousState() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { if (it->second->hasHistory()) it->second->saveCurrentValues(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ // void Material::updateResidual(GhostType ghost_type) { // AKANTU_DEBUG_IN(); // computeAllStresses(ghost_type); // assembleResidual(ghost_type); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); if (!finite_deformation) { auto & internal_force = const_cast &>(model.getInternalForce()); Mesh & mesh = fem.getMesh(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type); for (; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element) { const Array & shapes_derivatives = fem.getShapesDerivatives(*it, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = fem.getNbIntegrationPoints(*it, ghost_type); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by /// @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Array * sigma_dphi_dx = new Array(nb_element * nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Array * shapesd_filtered = new Array(0, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, *it, ghost_type, elem_filter); Array & stress_vect = this->stress(*it, ghost_type); Array::matrix_iterator sigma = stress_vect.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator B = shapesd_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_sigma_it = sigma_dphi_dx->begin(spatial_dimension, nb_nodes_per_element); for (UInt q = 0; q < nb_element * nb_quadrature_points; ++q, ++sigma, ++B, ++Bt_sigma_it) Bt_sigma_it->mul(*sigma, *B); delete shapesd_filtered; /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by * @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ Array * int_sigma_dphi_dx = new Array( nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); fem.integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, *it, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble model.getDOFManager().assembleElementalArrayLocalArray( *int_sigma_dphi_dx, internal_force, *it, ghost_type, -1, elem_filter); delete int_sigma_dphi_dx; } } } else { switch (spatial_dimension) { case 1: this->assembleInternalForces<1>(ghost_type); break; case 2: this->assembleInternalForces<2>(ghost_type); break; case 3: this->assembleInternalForces<3>(ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the gradu * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (const auto & type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, spatial_dimension, type, ghost_type, elem_filter); gradu_vect -= eigengradu(type, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllCauchyStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(finite_deformation, "The Cauchy stress can only be " "computed if you are working in " "finite deformation."); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { switch (spatial_dimension) { case 1: this->computeCauchyStress<1>(type, ghost_type); break; case 2: this->computeCauchyStress<2>(type, ghost_type); break; case 3: this->computeCauchyStress<3>(type, ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeCauchyStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::matrix_iterator gradu_it = this->gradu(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator gradu_end = this->gradu(el_type, ghost_type).end(dim, dim); Array::matrix_iterator piola_it = this->piola_kirchhoff_2(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator stress_it = this->stress(el_type, ghost_type).begin(dim, dim); Matrix F_tensor(dim, dim); for (; gradu_it != gradu_end; ++gradu_it, ++piola_it, ++stress_it) { Matrix & grad_u = *gradu_it; Matrix & piola = *piola_it; Matrix & sigma = *stress_it; gradUToF(grad_u, F_tensor); this->computeCauchyStressOnQuad(F_tensor, piola, sigma); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & displacement = model.getDisplacement(); // resizeInternalArray(gradu); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(displacement, gradu_vect, spatial_dimension, type, ghost_type, elem_filter); setToSteadyState(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { if (finite_deformation) { switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL<1>(type, ghost_type); assembleStiffnessMatrixL2<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrixNL<2>(type, ghost_type); assembleStiffnessMatrixL2<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrixNL<3>(type, ghost_type); assembleStiffnessMatrixL2<3>(type, ghost_type); break; } } } else { switch (spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(type, ghost_type); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size()) { const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array( nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapesd_filtered = new Array( nb_element, dim * nb_nodes_per_element, "filtered shapesd"); FEEngine::filterElementalData(fem.getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator Bt_D_B_it = bt_d_b->begin(dim * nb_nodes_per_element, dim * nb_nodes_per_element); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::matrix_iterator D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it) { Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_D.mul(B, D); Bt_D_B.mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); // Array & gradu_vect = delta_gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); // gradu_vect.resize(nb_quadrature_points * nb_element); // fem.gradientOnIntegrationPoints(model.getIncrement(), gradu_vect, // dim, type, ghost_type, &elem_filter); Array * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_s_b_size = dim * nb_nodes_per_element; Array * bt_s_b = new Array(nb_element * nb_quadrature_points, bt_s_b_size * bt_s_b_size, "B^t*D*B"); UInt piola_matrix_size = getCauchyStressMatrixSize(dim); Matrix B(piola_matrix_size, bt_s_b_size); Matrix Bt_S(bt_s_b_size, piola_matrix_size); Matrix S(piola_matrix_size, piola_matrix_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); Array::matrix_iterator Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); Array::matrix_iterator piola_it = piola_kirchhoff_2(type, ghost_type).begin(dim, dim); for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) { Matrix & Bt_S_B = *Bt_S_B_it; Matrix & Piola_kirchhoff_matrix = *piola_it; setCauchyStressMatrix(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.mul(B, S); Bt_S_B.mul(Bt_S, B); } delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_s_b_size * bt_s_b_size, "K_e"); fem.integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, elem_filter); delete bt_s_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_D_B_it = bt_d_b->begin(dim * nb_nodes_per_element, dim * nb_nodes_per_element); Array::matrix_iterator grad_u_it = gradu_vect.begin(dim, dim); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::matrix_iterator D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) { Matrix & grad_u = *grad_u_it; Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; // transferBMatrixToBL1 (*shapes_derivatives_filtered_it, B, // nb_nodes_per_element); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element); B += B2; Bt_D.mul(B, D); Bt_D_B.mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); Mesh & mesh = fem.getMesh(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array * shapesd_filtered = new Array( nb_element, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); // Set stress vectors UInt stress_size = getTangentStiffnessVoigtSize(dim); // Set matrices B and BNL* UInt bt_s_size = dim * nb_nodes_per_element; Array * bt_s = new Array(nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); Array::matrix_iterator grad_u_it = this->gradu(type, ghost_type).begin(dim, dim); Array::matrix_iterator grad_u_end = this->gradu(type, ghost_type).end(dim, dim); Array::matrix_iterator stress_it = this->piola_kirchhoff_2(type, ghost_type).begin(dim, dim); shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1); Matrix S_vect(stress_size, 1); Matrix B_tensor(stress_size, bt_s_size); Matrix B2_tensor(stress_size, bt_s_size); for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) { Matrix & grad_u = *grad_u_it; Matrix & r_it = *bt_s_it; Matrix & S_it = *stress_it; setCauchyStressArray(S_it, S_vect); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element); B_tensor += B2_tensor; r_it.mul(B_tensor, S_vect); } delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * r_e = new Array(nb_element, bt_s_size, "r_e"); fem.integrate(*bt_s, *r_e, bt_s_size, type, ghost_type, elem_filter); delete bt_s; model.getDOFManager().assembleElementalArrayLocalArray( *r_e, internal_force, type, ghost_type, -1., elem_filter); delete r_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { switch (spatial_dimension) { case 1: { computeAllStressesFromTangentModuli<1>(type, ghost_type); break; } case 2: { computeAllStressesFromTangentModuli<2>(type, ghost_type); break; } case 3: { computeAllStressesFromTangentModuli<3>(type, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); Array & disp = model.getDisplacement(); fem.gradientOnIntegrationPoints(disp, gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_moduli_size = getTangentStiffnessVoigtSize(dim); Array * tangent_moduli_tensors = new Array( nb_element * nb_quadrature_points, tangent_moduli_size * tangent_moduli_size, "tangent_moduli_tensors"); tangent_moduli_tensors->clear(); computeTangentModuli(type, *tangent_moduli_tensors, ghost_type); Array * shapesd_filtered = new Array( nb_element, dim * nb_nodes_per_element, "filtered shapesd"); FEEngine::filterElementalData(fem.getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array filtered_u(nb_element, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(fem.getMesh(), disp, filtered_u, type, ghost_type, elem_filter); /// compute @f$\mathbf{D} \mathbf{B} \mathbf{u}@f$ Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator D_it = tangent_moduli_tensors->begin(tangent_moduli_size, tangent_moduli_size); Array::matrix_iterator sigma_it = stress(type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::vector_iterator u_it = filtered_u.begin(spatial_dimension * nb_nodes_per_element); Matrix B(tangent_moduli_size, spatial_dimension * nb_nodes_per_element); Vector Bu(tangent_moduli_size); Vector DBu(tangent_moduli_size); for (UInt e = 0; e < nb_element; ++e, ++u_it) { for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it, ++shapes_derivatives_filtered_it, ++sigma_it) { Vector & u = *u_it; Matrix & sigma = *sigma_it; Matrix & D = *D_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bu.mul(B, u); DBu.mul(D, Bu); // Voigt notation to full symmetric tensor for (UInt i = 0; i < dim; ++i) sigma(i, i) = DBu(i); if (dim == 2) { sigma(0, 1) = sigma(1, 0) = DBu(2); } else if (dim == 3) { sigma(1, 2) = sigma(2, 1) = DBu(3); sigma(0, 2) = sigma(2, 0) = DBu(4); sigma(0, 1) = sigma(1, 0) = DBu(5); } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElements() { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { computePotentialEnergy(type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType, GhostType) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElements(); /// integrate the potential energy for each type of elements for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { epot += fem.integrate(potential_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy(ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real epot = 0.; Vector epot_on_quad_points(fem.getNbIntegrationPoints(type)); computePotentialEnergyByElement(type, index, epot_on_quad_points); epot = fem.integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ -Real Material::getEnergy(std::string type) { +Real Material::getEnergy(const std::string & type) { AKANTU_DEBUG_IN(); if (type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ -Real Material::getEnergy(std::string energy_id, ElementType type, UInt index) { +Real Material::getEnergy(const std::string & energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "potential") return getPotentialEnergy(type, index); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation( const ElementTypeMapArray & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); this->fem.initElementalFieldInterpolationFromIntegrationPoints( interpolation_points_coordinates, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, &(this->element_filter)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type) { this->fem.interpolateElementalFieldFromIntegrationPoints( this->stress, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, result, ghost_type, &(this->element_filter)); } /* -------------------------------------------------------------------------- */ void Material::interpolateStressOnFacets( ElementTypeMapArray & result, ElementTypeMapArray & by_elem_result, const GhostType ghost_type) { interpolateStress(by_elem_result, ghost_type); UInt stress_size = this->stress.getNbComponent(); const Mesh & mesh = this->model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_fil = element_filter(type, ghost_type); Array & by_elem_res = by_elem_result(type, ghost_type); UInt nb_element = elem_fil.size(); UInt nb_element_full = this->model.getMesh().getNbElement(type, ghost_type); UInt nb_interpolation_points_per_elem = by_elem_res.size() / nb_element_full; const Array & facet_to_element = mesh_facets.getSubelementToElement(type, ghost_type); ElementType type_facet = Mesh::getFacetType(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); UInt nb_quad_per_facet = nb_interpolation_points_per_elem / nb_facet_per_elem; Element element_for_comparison{type, 0, ghost_type}; const Array> * element_to_facet = nullptr; GhostType current_ghost_type = _casper; Array * result_vec = nullptr; Array::const_matrix_iterator result_it = by_elem_res.begin_reinterpret( stress_size, nb_interpolation_points_per_elem, nb_element_full); for (UInt el = 0; el < nb_element; ++el) { UInt global_el = elem_fil(el); element_for_comparison.element = global_el; for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element facet_elem = facet_to_element(global_el, f); UInt global_facet = facet_elem.element; if (facet_elem.ghost_type != current_ghost_type) { current_ghost_type = facet_elem.ghost_type; element_to_facet = &mesh_facets.getElementToSubelement( type_facet, current_ghost_type); result_vec = &result(type_facet, current_ghost_type); } bool is_second_element = (*element_to_facet)(global_facet)[0] != element_for_comparison; for (UInt q = 0; q < nb_quad_per_facet; ++q) { Vector result_local(result_vec->storage() + (global_facet * nb_quad_per_facet + q) * result_vec->getNbComponent() + is_second_element * stress_size, stress_size); const Matrix & result_tmp(result_it[global_el]); result_local = result_tmp(f * nb_quad_per_facet + q); } } } } } /* -------------------------------------------------------------------------- */ template const Array & Material::getArray(__attribute__((unused)) const ID & vect_id, __attribute__((unused)) const ElementType & type, __attribute__((unused)) const GhostType & ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template Array & Material::getArray(__attribute__((unused)) const ID & vect_id, __attribute__((unused)) const ElementType & type, __attribute__((unused)) const GhostType & ghost_type) { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template const InternalField & Material::getInternal(__attribute__((unused)) const ID & int_id) const { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template InternalField & Material::getInternal(__attribute__((unused)) const ID & int_id) { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { auto it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { auto it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ void Material::addElements(const Array & elements_to_add) { AKANTU_DEBUG_IN(); UInt mat_id = model.getInternalIndexFromID(getID()); Array::const_iterator el_begin = elements_to_add.begin(); Array::const_iterator el_end = elements_to_add.end(); for (; el_begin != el_end; ++el_begin) { const Element & element = *el_begin; Array & mat_indexes = model.getMaterialByElement(element.type, element.ghost_type); Array & mat_loc_num = model.getMaterialLocalNumbering(element.type, element.ghost_type); UInt index = this->addElement(element.type, element.element, element.ghost_type); mat_indexes(element.element) = mat_id; mat_loc_num(element.element) = index; } this->resizeInternals(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::removeElements(const Array & elements_to_remove) { AKANTU_DEBUG_IN(); Array::const_iterator el_begin = elements_to_remove.begin(); Array::const_iterator el_end = elements_to_remove.end(); if (el_begin == el_end) return; ElementTypeMapArray material_local_new_numbering( "remove mat filter elem", getID(), getMemoryID()); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; element.ghost_type = ghost_type; ElementTypeMapArray::type_iterator it = element_filter.firstType(_all_dimensions, ghost_type, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, ghost_type, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; element.type = type; Array & elem_filter = this->element_filter(type, ghost_type); Array & mat_loc_num = this->model.getMaterialLocalNumbering(type, ghost_type); if (!material_local_new_numbering.exists(type, ghost_type)) material_local_new_numbering.alloc(elem_filter.size(), 1, type, ghost_type); Array & mat_renumbering = material_local_new_numbering(type, ghost_type); UInt nb_element = elem_filter.size(); Array elem_filter_tmp; UInt new_id = 0; for (UInt el = 0; el < nb_element; ++el) { element.element = elem_filter(el); if (std::find(el_begin, el_end, element) == el_end) { elem_filter_tmp.push_back(element.element); mat_renumbering(el) = new_id; mat_loc_num(element.element) = new_id; ++new_id; } else { mat_renumbering(el) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.size()); elem_filter.copy(elem_filter_tmp); } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resizeInternals() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::onElementsAdded(const Array &, const NewElementsEvent &) { this->resizeInternals(); } /* -------------------------------------------------------------------------- */ void Material::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { UInt my_num = model.getInternalIndexFromID(getID()); ElementTypeMapArray material_local_new_numbering( "remove mat filter elem", getID(), getMemoryID()); Array::const_iterator el_begin = element_list.begin(); Array::const_iterator el_end = element_list.end(); for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType gt = *g; ElementTypeMapArray::type_iterator it = new_numbering.firstType(_all_dimensions, gt, _ek_not_defined); ElementTypeMapArray::type_iterator end = new_numbering.lastType(_all_dimensions, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if (element_filter.exists(type, gt) && element_filter(type, gt).size()) { Array & elem_filter = element_filter(type, gt); Array & mat_indexes = this->model.getMaterialByElement(*it, gt); Array & mat_loc_num = this->model.getMaterialLocalNumbering(*it, gt); UInt nb_element = this->model.getMesh().getNbElement(type, gt); // all materials will resize of the same size... mat_indexes.resize(nb_element); mat_loc_num.resize(nb_element); if (!material_local_new_numbering.exists(type, gt)) material_local_new_numbering.alloc(elem_filter.size(), 1, type, gt); Array & mat_renumbering = material_local_new_numbering(type, gt); const Array & renumbering = new_numbering(type, gt); Array elem_filter_tmp; UInt ni = 0; Element el{type, 0, gt}; for (UInt i = 0; i < elem_filter.size(); ++i) { el.element = elem_filter(i); if (std::find(el_begin, el_end, el) == el_end) { UInt new_el = renumbering(el.element); AKANTU_DEBUG_ASSERT( new_el != UInt(-1), "A not removed element as been badly renumbered"); elem_filter_tmp.push_back(new_el); mat_renumbering(i) = ni; mat_indexes(new_el) = my_num; mat_loc_num(new_el) = ni; ++ni; } else { mat_renumbering(i) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.size()); elem_filter.copy(elem_filter_tmp); } } } for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); } /* -------------------------------------------------------------------------- */ void Material::beforeSolveStep() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::afterSolveStep() { for (auto & type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { this->updateEnergies(type, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDamageIteration() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onDamageUpdate() { ElementTypeMapArray::type_iterator it = this->element_filter.firstType( _all_dimensions, _not_ghost, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined); for (; it != end; ++it) { this->updateEnergiesAfterDamage(*it, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDump() { if (this->isFiniteDeformation()) this->computeAllCauchyStresses(_not_ghost); } /* -------------------------------------------------------------------------- */ void Material::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /// extrapolate internal values void Material::extrapolateInternal(const ID & id, const Element & element, __attribute__((unused)) const Matrix & point, Matrix & extrapolated) { if (this->isInternal(id, element.kind())) { UInt nb_element = this->element_filter(element.type, element.ghost_type).size(); const ID name = this->getID() + ":" + id; UInt nb_quads = this->internal_vectors_real[name]->getFEEngine().getNbIntegrationPoints( element.type, element.ghost_type); const Array & internal = this->getArray(id, element.type, element.ghost_type); UInt nb_component = internal.getNbComponent(); Array::const_matrix_iterator internal_it = internal.begin_reinterpret(nb_component, nb_quads, nb_element); Element local_element = this->convertToLocalElement(element); /// instead of really extrapolating, here the value of the first GP /// is copied into the result vector. This works only for linear /// elements /// @todo extrapolate!!!! AKANTU_DEBUG_WARNING("This is a fix, values are not truly extrapolated"); const Matrix & values = internal_it[local_element.element]; UInt index = 0; Vector tmp(nb_component); for (UInt j = 0; j < values.cols(); ++j) { tmp = values(j); if (tmp.norm() > 0) { index = j; break; } } for (UInt i = 0; i < extrapolated.size(); ++i) { extrapolated(i) = values(index); } } else { Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); extrapolated = default_values; } } /* -------------------------------------------------------------------------- */ void Material::applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const GhostType ghost_type) { for (auto && type : element_filter.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) { if (!element_filter(type, ghost_type).size()) continue; auto eigen_it = this->eigengradu(type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto eigen_end = this->eigengradu(type, ghost_type) .end(spatial_dimension, spatial_dimension); for (; eigen_it != eigen_end; ++eigen_it) { auto & current_eigengradu = *eigen_it; current_eigengradu = prescribed_eigen_grad_u; } } } } // namespace akantu diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index db578d06c..a99d4cb27 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,681 +1,681 @@ /** * @file material.hh * * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Nov 25 2015 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_factory.hh" #include "aka_voigthelper.hh" #include "data_accessor.hh" #include "integration_point.hh" #include "parsable.hh" #include "parser.hh" /* -------------------------------------------------------------------------- */ #include "internal_field.hh" #include "random_internal_field.hh" /* -------------------------------------------------------------------------- */ #include "mesh_events.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; } // namespace akantu namespace akantu { /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = * ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(const ElementType & el_type, * Array & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : public Memory, public DataAccessor, public Parsable, public MeshEventHandler, protected SolidMechanicsModelEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: #if __cplusplus > 199711L Material(const Material & mat) = delete; Material & operator=(const Material & mat) = delete; #endif /// Initialize material with defaults Material(SolidMechanicsModel & model, const ID & id = ""); /// Initialize material with custom mesh & fe_engine Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); /// Destructor ~Material() override; protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Function that materials can/should reimplement */ /* ------------------------------------------------------------------------ */ protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the tangent stiffness matrix virtual void computeTangentModuli(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the potential energy for an element virtual void computePotentialEnergyByElement(__attribute__((unused)) ElementType type, __attribute__((unused)) UInt index, __attribute__((unused)) Vector & epot_on_quad_points) { AKANTU_DEBUG_TO_IMPLEMENT(); } virtual void updateEnergies(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} virtual void updateEnergiesAfterDamage(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} /// set the material to steady state (to be implemented for materials that /// need it) virtual void setToSteadyState(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} /// function called to update the internal parameters when the modifiable /// parameters are modified virtual void updateInternalParameters() {} public: /// extrapolate internal values virtual void extrapolateInternal(const ID & id, const Element & element, const Matrix & points, Matrix & extrapolated); /// compute the p-wave speed in the material virtual Real getPushWaveSpeed(__attribute__((unused)) const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the s-wave speed in the material virtual Real getShearWaveSpeed(__attribute__((unused)) const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// get a material celerity to compute the stable time step (default: is the /// push wave speed) virtual Real getCelerity(const Element & element) const { return getPushWaveSpeed(element); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void registerInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } template void unregisterInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material // virtual void updateResidual(GhostType ghost_type = _not_ghost); /// assemble the residual for this material virtual void assembleInternalForces(GhostType ghost_type); /// save the stress in the previous_stress if needed virtual void savePreviousState(); /// compute the stresses for this material virtual void computeAllStresses(GhostType ghost_type = _not_ghost); virtual void computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost); virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); inline UInt addElement(const Element & element); /// add many elements at once void addElements(const Array & elements_to_add); /// remove many element at once void removeElements(const Array & elements_to_remove); /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points */ void interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type = _not_ghost); /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points and store the * results per facet */ void interpolateStressOnFacets(ElementTypeMapArray & result, ElementTypeMapArray & by_elem_result, const GhostType ghost_type = _not_ghost); /** * function to initialize the elemental field interpolation * function by inverting the quadrature points' coordinates */ void initElementalFieldInterpolation( const ElementTypeMapArray & interpolation_points_coordinates); /* ------------------------------------------------------------------------ */ /* Common part */ /* ------------------------------------------------------------------------ */ protected: /* ------------------------------------------------------------------------ */ inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// compute the potential energy by element void computePotentialEnergyByElements(); /// resize the intenals arrays virtual void resizeInternals(); /* ------------------------------------------------------------------------ */ /* Finite deformation functions */ /* This functions area implementing what is described in the paper of Bathe */ /* et al, in IJNME, Finite Element Formulations For Large Deformation */ /* Dynamic Analysis, Vol 9, 353-386, 1975 */ /* ------------------------------------------------------------------------ */ protected: /// assemble the residual template void assembleInternalForces(GhostType ghost_type); /// Computation of Cauchy stress tensor in the case of finite deformation from /// the 2nd Piola-Kirchhoff for a given element type template void computeCauchyStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation /// gradient template inline void computeCauchyStressOnQuad(const Matrix & F, const Matrix & S, Matrix & cauchy, const Real & C33 = 1.0) const; template void computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// assembling in finite deformation template void assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type); /// Size of the Stress matrix for the case of finite deformation see: Bathe et /// al, IJNME, Vol 9, 353-386, 1975 inline UInt getCauchyStressMatrixSize(UInt spatial_dimension) const; /// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386, /// 1975 template inline void setCauchyStressMatrix(const Matrix & S_t, Matrix & sigma); /// write the stress tensor in the Voigt notation. template inline void setCauchyStressArray(const Matrix & S_t, Matrix & sigma_voight); /* ------------------------------------------------------------------------ */ /* Conversion functions */ /* ------------------------------------------------------------------------ */ public: template static inline void gradUToF(const Matrix & grad_u, Matrix & F); static inline void rightCauchy(const Matrix & F, Matrix & C); static inline void leftCauchy(const Matrix & F, Matrix & B); template static inline void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); template static inline void gradUToGreenStrain(const Matrix & grad_u, Matrix & epsilon); static inline Real stressToVonMises(const Matrix & stress); protected: /// converts global element to local element inline Element convertToLocalElement(const Element & global_element) const; /// converts local element to global element inline Element convertToGlobalElement(const Element & local_element) const; /// converts global quadrature point to local quadrature point inline IntegrationPoint convertToLocalPoint(const IntegrationPoint & global_point) const; /// converts local quadrature point to global quadrature point inline IntegrationPoint convertToGlobalPoint(const IntegrationPoint & local_point) const; /* ------------------------------------------------------------------------ */ /* DataAccessor 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; template inline void packElementDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); /* ------------------------------------------------------------------------ */ /* MeshEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ void onNodesAdded(const Array &, const NewNodesEvent &) override{}; void onNodesRemoved(const Array &, const Array &, const RemovedNodesEvent &) override{}; void onElementsAdded(const Array & element_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{}; /* ------------------------------------------------------------------------ */ /* SolidMechanicsModelEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: virtual void beforeSolveStep(); virtual void afterSolveStep(); void onDamageIteration() override; void onDamageUpdate() override; void onDump() override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Name, name, const std::string &); AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, Memory::getID(), const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// return the potential energy for the subset of elements contained by the /// material Real getPotentialEnergy(); /// return the potential energy for the provided element Real getPotentialEnergy(ElementType & type, UInt index); /// return the energy (identified by id) for the subset of elements contained /// by the material - virtual Real getEnergy(std::string energy_id); + virtual Real getEnergy(const std::string & energy_id); /// return the energy (identified by id) for the provided element - virtual Real getEnergy(std::string energy_id, ElementType type, UInt index); + virtual Real getEnergy(const std::string & energy_id, ElementType type, UInt index); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray &); AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray &); AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); AKANTU_GET_MACRO(FEEngine, fem, FEEngine &); bool isNonLocal() const { return is_non_local; } template const Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; template Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); template const InternalField & getInternal(const ID & id) const; template InternalField & getInternal(const ID & id); template inline bool isInternal(const ID & id, const ElementKind & element_kind) const; template ElementTypeMap getInternalDataPerElem(const ID & id, const ElementKind & element_kind) const; bool isFiniteDeformation() const { return finite_deformation; } bool isInelasticDeformation() const { return inelastic_deformation; } template inline void setParam(const ID & param, T value); inline const Parameter & getParam(const ID & param) const; template void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) const; /// apply a constant eigengrad_u everywhere in the material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const GhostType = _not_ghost); /// specify if the matrix need to be recomputed for this material virtual bool hasStiffnessMatrixChanged() { return true; } protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// boolean to know if the material has been initialized bool is_init; std::map *> internal_vectors_real; std::map *> internal_vectors_uint; std::map *> internal_vectors_bool; protected: /// Link to the fem object in the model FEEngine & fem; /// Finite deformation bool finite_deformation; /// Finite deformation bool inelastic_deformation; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel & model; /// density : rho Real rho; /// spatial dimension UInt spatial_dimension; /// list of element handled by the material ElementTypeMapArray element_filter; /// stresses arrays ordered by element types InternalField stress; /// eigengrad_u arrays ordered by element types InternalField eigengradu; /// grad_u arrays ordered by element types InternalField gradu; /// Green Lagrange strain (Finite deformation) InternalField green_strain; /// Second Piola-Kirchhoff stress tensor arrays ordered by element types /// (Finite deformation) InternalField piola_kirchhoff_2; /// potential energy by element InternalField potential_energy; /// tell if using in non local mode or not bool is_non_local; /// tell if the material need the previous stress state bool use_previous_stress; /// tell if the material need the previous strain state bool use_previous_gradu; /// elemental field interpolation coordinates InternalField interpolation_inverse_coordinates; /// elemental field interpolation points InternalField interpolation_points_matrices; /// vector that contains the names of all the internals that need to /// be transferred when material interfaces move std::vector internals_to_transfer; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "material_inline_impl.cc" #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ /// This can be used to automatically write the loop on quadrature points in /// functions such as computeStress. This macro in addition to write the loop /// provides two tensors (matrices) sigma and grad_u #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ Array::matrix_iterator gradu_it = \ this->gradu(el_type, ghost_type) \ .begin(this->spatial_dimension, this->spatial_dimension); \ Array::matrix_iterator gradu_end = \ this->gradu(el_type, ghost_type) \ .end(this->spatial_dimension, this->spatial_dimension); \ \ this->stress(el_type, ghost_type) \ .resize(this->gradu(el_type, ghost_type).size()); \ \ Array::iterator> stress_it = \ this->stress(el_type, ghost_type) \ .begin(this->spatial_dimension, this->spatial_dimension); \ \ if (this->isFiniteDeformation()) { \ this->piola_kirchhoff_2(el_type, ghost_type) \ .resize(this->gradu(el_type, ghost_type).size()); \ stress_it = this->piola_kirchhoff_2(el_type, ghost_type) \ .begin(this->spatial_dimension, this->spatial_dimension); \ } \ \ for (; gradu_it != gradu_end; ++gradu_it, ++stress_it) { \ Matrix & __attribute__((unused)) grad_u = *gradu_it; \ Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END } /// This can be used to automatically write the loop on quadrature points in /// functions such as computeTangentModuli. This macro in addition to write the /// loop provides two tensors (matrices) sigma_tensor, grad_u, and a matrix /// where the elemental tangent moduli should be stored in Voigt Notation #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Array::matrix_iterator gradu_it = \ this->gradu(el_type, ghost_type) \ .begin(this->spatial_dimension, this->spatial_dimension); \ Array::matrix_iterator gradu_end = \ this->gradu(el_type, ghost_type) \ .end(this->spatial_dimension, this->spatial_dimension); \ Array::matrix_iterator sigma_it = \ this->stress(el_type, ghost_type) \ .begin(this->spatial_dimension, this->spatial_dimension); \ \ tangent_mat.resize(this->gradu(el_type, ghost_type).size()); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(this->spatial_dimension); \ Array::matrix_iterator tangent_it = \ tangent_mat.begin(tangent_size, tangent_size); \ \ for (; gradu_it != gradu_end; ++gradu_it, ++sigma_it, ++tangent_it) { \ Matrix & __attribute__((unused)) grad_u = *gradu_it; \ Matrix & __attribute__((unused)) sigma_tensor = *sigma_it; \ Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END } /* -------------------------------------------------------------------------- */ namespace akantu { using MaterialFactory = Factory; } // namespace akantu #define INSTANTIATE_MATERIAL_ONLY(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #define MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name) \ [](UInt dim, const ID &, SolidMechanicsModel & model, \ const ID & id) -> std::unique_ptr { \ switch (dim) { \ case 1: \ return std::make_unique>(model, id); \ case 2: \ return std::make_unique>(model, id); \ case 3: \ return std::make_unique>(model, id); \ default: \ AKANTU_EXCEPTION("The dimension " \ << dim << "is not a valid dimension for the material " \ << #id); \ } \ } #define INSTANTIATE_MATERIAL(id, mat_name) \ INSTANTIATE_MATERIAL_ONLY(mat_name); \ static bool material_is_alocated_##id = \ MaterialFactory::getInstance().registerAllocator( \ #id, MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name)) #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_damage.hh b/src/model/solid_mechanics/materials/material_damage/material_damage.hh index 475c6d0b5..fb72e1549 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_damage.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_damage.hh @@ -1,112 +1,113 @@ /** * @file material_damage.hh * * @author Marion Estelle Chambart * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Sat Jul 11 2015 * * @brief Material isotropic elastic * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DAMAGE_HH__ #define __AKANTU_MATERIAL_DAMAGE_HH__ namespace akantu { template class Parent = MaterialElastic> class MaterialDamage : public Parent { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialDamage(SolidMechanicsModel & model, const ID & id = ""); ~MaterialDamage() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial() override; /// compute the tangent stiffness matrix for an element type void computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override; bool hasStiffnessMatrixChanged() override { return true; } protected: /// update the dissipated energy, must be called after the stress have been /// computed void updateEnergies(ElementType el_type, GhostType ghost_type) override; /// compute the tangent stiffness matrix for a given quadrature point inline void computeTangentModuliOnQuad(Matrix & tangent, Real & dam); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// give the dissipated energy for the time step Real getDissipatedEnergy() const; - Real getEnergy(std::string type) override; - Real getEnergy(std::string energy_id, ElementType type, UInt index) override { + Real getEnergy(const std::string & type) override; + Real getEnergy(const std::string & energy_id, ElementType type, + UInt index) override { return Parent::getEnergy(energy_id, type, index); }; AKANTU_GET_MACRO_NOT_CONST(Damage, damage, ElementTypeMapArray &); AKANTU_GET_MACRO(Damage, damage, const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real) /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// damage internal variable InternalField damage; /// dissipated energy InternalField dissipated_energy; /// contain the current value of @f$ \int_0^{\epsilon}\sigma(\omega)d\omega /// @f$ the dissipated energy InternalField int_sigma; }; -} // akantu +} // namespace akantu #include "material_damage_tmpl.hh" #endif /* __AKANTU_MATERIAL_DAMAGE_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh b/src/model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh index bbe847c92..2876eb704 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh @@ -1,180 +1,181 @@ /** * @file material_damage_tmpl.hh * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Aug 18 2015 * * @brief Specialization of the material class for the damage material * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_damage.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class Parent> MaterialDamage::MaterialDamage( SolidMechanicsModel & model, const ID & id) : Parent(model, id), damage("damage", *this), dissipated_energy("damage dissipated energy", *this), int_sigma("integral of sigma", *this) { AKANTU_DEBUG_IN(); this->is_non_local = false; this->use_previous_stress = true; this->use_previous_gradu = true; this->damage.initialize(1); this->dissipated_energy.initialize(1); this->int_sigma.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialDamage::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the dissipated energy in each element by a trapezoidal approximation * of * @f$ Ed = \int_0^{\epsilon}\sigma(\omega)d\omega - * \frac{1}{2}\sigma:\epsilon@f$ */ template class Parent> void MaterialDamage::updateEnergies( ElementType el_type, GhostType ghost_type) { Parent::updateEnergies(el_type, ghost_type); this->computePotentialEnergy(el_type, ghost_type); Array::matrix_iterator epsilon_p = this->gradu.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); Array::matrix_iterator sigma_p = this->stress.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); Array::const_scalar_iterator epot = this->potential_energy(el_type, ghost_type).begin(); Array::scalar_iterator ints = this->int_sigma(el_type, ghost_type).begin(); Array::scalar_iterator ed = this->dissipated_energy(el_type, ghost_type).begin(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix delta_gradu_it(*gradu_it); delta_gradu_it -= *epsilon_p; Matrix sigma_h(sigma); sigma_h += *sigma_p; Real dint = .5 * sigma_h.doubleDot(delta_gradu_it); *ints += dint; *ed = *ints - *epot; ++epsilon_p; ++sigma_p; ++epot; ++ints; ++ed; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialDamage::computeTangentModuli( const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); Parent::computeTangentModuli(el_type, tangent_matrix, ghost_type); Real * dam = this->damage(el_type, ghost_type).storage(); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); computeTangentModuliOnQuad(tangent, *dam); ++dam; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialDamage::computeTangentModuliOnQuad( Matrix & tangent, Real & dam) { tangent *= (1 - dam); } /* -------------------------------------------------------------------------- */ template class Parent> Real MaterialDamage::getDissipatedEnergy() const { AKANTU_DEBUG_IN(); Real de = 0.; /// integrate the dissipated energy for each type of elements for (auto & type : this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { de += this->fem.integrate(dissipated_energy(type, _not_ghost), type, _not_ghost, this->element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return de; } /* -------------------------------------------------------------------------- */ template class Parent> -Real MaterialDamage::getEnergy(std::string type) { +Real MaterialDamage::getEnergy( + const std::string & type) { if (type == "dissipated") return getDissipatedEnergy(); else return Parent::getEnergy(type); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_plastic/material_plastic.cc b/src/model/solid_mechanics/materials/material_plastic/material_plastic.cc index 2f407b6e0..bfe07f012 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_plastic.cc +++ b/src/model/solid_mechanics/materials/material_plastic/material_plastic.cc @@ -1,209 +1,209 @@ /** * @file material_plastic.cc * * @author Lucas Frerot * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Mon Apr 07 2014 * @date last modification: Tue Aug 18 2015 * * @brief Implemantation of the akantu::MaterialPlastic class * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_plastic.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialPlastic::MaterialPlastic(SolidMechanicsModel & model, const ID & id) : MaterialElastic(model, id), iso_hardening("iso_hardening", *this), inelastic_strain("inelastic_strain", *this), plastic_energy("plastic_energy", *this), d_plastic_energy("d_plastic_energy", *this) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } template MaterialPlastic::MaterialPlastic(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : MaterialElastic(model, dim, mesh, fe_engine, id), iso_hardening("iso_hardening", *this, dim, fe_engine, this->element_filter), inelastic_strain("inelastic_strain", *this, dim, fe_engine, this->element_filter), plastic_energy("plastic_energy", *this, dim, fe_engine, this->element_filter), d_plastic_energy("d_plastic_energy", *this, dim, fe_engine, this->element_filter) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialPlastic::initialize() { this->registerParam("h", h, Real(0.), _pat_parsable | _pat_modifiable, "Hardening modulus"); this->registerParam("sigma_y", sigma_y, Real(0.), _pat_parsable | _pat_modifiable, "Yield stress"); this->iso_hardening.initialize(1); this->iso_hardening.initializeHistory(); this->plastic_energy.initialize(1); this->d_plastic_energy.initialize(1); this->use_previous_stress = true; this->use_previous_gradu = true; this->use_previous_stress_thermal = true; this->inelastic_strain.initialize(spatial_dimension * spatial_dimension); this->inelastic_strain.initializeHistory(); } /* -------------------------------------------------------------------------- */ template -Real MaterialPlastic::getEnergy(std::string type) { +Real MaterialPlastic::getEnergy(const std::string & type) { if (type == "plastic") return getPlasticEnergy(); else return MaterialElastic::getEnergy(type); return 0.; } /* -------------------------------------------------------------------------- */ template Real MaterialPlastic::getPlasticEnergy() { AKANTU_DEBUG_IN(); Real penergy = 0.; for (auto & type : this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { penergy += this->fem.integrate(plastic_energy(type, _not_ghost), type, _not_ghost, this->element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return penergy; } /* -------------------------------------------------------------------------- */ template void MaterialPlastic::computePotentialEnergy( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (ghost_type != _not_ghost) return; Array::scalar_iterator epot = this->potential_energy(el_type, ghost_type).begin(); Array::const_iterator> inelastic_strain_it = this->inelastic_strain(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix elastic_strain(spatial_dimension, spatial_dimension); elastic_strain.copy(grad_u); elastic_strain -= *inelastic_strain_it; MaterialElastic::computePotentialEnergyOnQuad( elastic_strain, sigma, *epot); ++epot; ++inelastic_strain_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialPlastic::updateEnergies(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialElastic::updateEnergies(el_type, ghost_type); Array::iterator<> pe_it = this->plastic_energy(el_type, ghost_type).begin(); Array::iterator<> wp_it = this->d_plastic_energy(el_type, ghost_type).begin(); Array::iterator> inelastic_strain_it = this->inelastic_strain(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); Array::iterator> previous_inelastic_strain_it = this->inelastic_strain.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); Array::matrix_iterator previous_sigma = this->stress.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix delta_strain_it(*inelastic_strain_it); delta_strain_it -= *previous_inelastic_strain_it; Matrix sigma_h(sigma); sigma_h += *previous_sigma; *wp_it = .5 * sigma_h.doubleDot(delta_strain_it); *pe_it += *wp_it; ++pe_it; ++wp_it; ++inelastic_strain_it; ++previous_inelastic_strain_it; ++previous_sigma; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL_ONLY(MaterialPlastic); } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_plastic/material_plastic.hh b/src/model/solid_mechanics/materials/material_plastic/material_plastic.hh index 2cff5aa8d..477b993b6 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_plastic.hh +++ b/src/model/solid_mechanics/materials/material_plastic/material_plastic.hh @@ -1,143 +1,143 @@ /** * @file material_plastic.hh * * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Oct 08 2015 * * @brief Common interface for plastic materials * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_PLASTIC_HH__ #define __AKANTU_MATERIAL_PLASTIC_HH__ namespace akantu { /** * Parent class for the plastic constitutive laws * parameters in the material files : * - h : Hardening parameter (default: 0) * - sigmay : Yield stress */ template class MaterialPlastic : public MaterialElastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialPlastic(SolidMechanicsModel & model, const ID & id = ""); MaterialPlastic(SolidMechanicsModel & model, UInt a_dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get the energy specifying the type for the time step - Real getEnergy(std::string type) override; + Real getEnergy(const std::string & type) override; /// Compute the plastic energy void updateEnergies(ElementType el_type, GhostType ghost_type = _not_ghost) override; /// Compute the true potential energy void computePotentialEnergy(ElementType el_type, GhostType ghost_type) override; protected: /// compute the stress and inelastic strain for the quadrature point inline void computeStressAndInelasticStrainOnQuad(const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelas_strain, const Matrix & previous_inelas_strain, const Matrix & delta_inelastic_strain) const; inline void computeStressAndInelasticStrainOnQuad(const Matrix & delta_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelas_strain, const Matrix & previous_inelas_strain, const Matrix & delta_inelastic_strain) const; /// get the plastic energy for the time step Real getPlasticEnergy(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Yield stresss Real sigma_y; /// hardening modulus Real h; /// isotropic hardening, r InternalField iso_hardening; /// inelastic strain arrays ordered by element types (inelastic deformation) InternalField inelastic_strain; /// Plastic energy InternalField plastic_energy; /// @todo : add a coefficient beta that will multiply the plastic energy increment /// to compute the energy converted to heat /// Plastic energy increment InternalField d_plastic_energy; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ } // akantu #include "material_plastic_inline_impl.cc" #endif /* __AKANTU_MATERIAL_PLASTIC_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc b/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc index 9aac29ef4..2ae5024da 100644 --- a/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc +++ b/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc @@ -1,295 +1,326 @@ /** * @file material_standard_linear_solid_deviatoric.cc * * @author David Simon Kammer * @author Nicolas Richart * @author Vladislav Yastrebov * * @date creation: Wed May 04 2011 * @date last modification: Thu Oct 15 2015 * * @brief Material Visco-elastic * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_standard_linear_solid_deviatoric.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -template -MaterialStandardLinearSolidDeviatoric::MaterialStandardLinearSolidDeviatoric(SolidMechanicsModel & model, const ID & id) : - MaterialElastic(model, id), - stress_dev("stress_dev", *this), - history_integral("history_integral", *this), - dissipated_energy("dissipated_energy", *this) { +template +MaterialStandardLinearSolidDeviatoric:: + MaterialStandardLinearSolidDeviatoric(SolidMechanicsModel & model, + const ID & id) + : MaterialElastic(model, id), + stress_dev("stress_dev", *this), + history_integral("history_integral", *this), + dissipated_energy("dissipated_energy", *this) { AKANTU_DEBUG_IN(); - this->registerParam("Eta", eta, Real(1.), _pat_parsable | _pat_modifiable, "Viscosity"); - this->registerParam("Ev", Ev, Real(1.), _pat_parsable | _pat_modifiable, "Stiffness of the viscous element"); - this->registerParam("Einf", E_inf, Real(1.), _pat_readable, "Stiffness of the elastic element"); + this->registerParam("Eta", eta, Real(1.), _pat_parsable | _pat_modifiable, + "Viscosity"); + this->registerParam("Ev", Ev, Real(1.), _pat_parsable | _pat_modifiable, + "Stiffness of the viscous element"); + this->registerParam("Einf", E_inf, Real(1.), _pat_readable, + "Stiffness of the elastic element"); UInt stress_size = spatial_dimension * spatial_dimension; - this->stress_dev .initialize(stress_size); - this->history_integral .initialize(stress_size); + this->stress_dev.initialize(stress_size); + this->history_integral.initialize(stress_size); this->dissipated_energy.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template +template void MaterialStandardLinearSolidDeviatoric::initMaterial() { AKANTU_DEBUG_IN(); updateInternalParameters(); MaterialElastic::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialStandardLinearSolidDeviatoric::updateInternalParameters() { +template +void MaterialStandardLinearSolidDeviatoric< + spatial_dimension>::updateInternalParameters() { MaterialElastic::updateInternalParameters(); E_inf = this->E - this->Ev; } /* -------------------------------------------------------------------------- */ -template -void MaterialStandardLinearSolidDeviatoric::setToSteadyState(ElementType el_type, GhostType ghost_type) { +template +void MaterialStandardLinearSolidDeviatoric::setToSteadyState( + ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - Array & stress_dev_vect = stress_dev(el_type, ghost_type); + Array & stress_dev_vect = stress_dev(el_type, ghost_type); Array & history_int_vect = history_integral(el_type, ghost_type); - Array::matrix_iterator stress_d = stress_dev_vect.begin(spatial_dimension, spatial_dimension); - Array::matrix_iterator history_int = history_int_vect.begin(spatial_dimension, spatial_dimension); + Array::matrix_iterator stress_d = + stress_dev_vect.begin(spatial_dimension, spatial_dimension); + Array::matrix_iterator history_int = + history_int_vect.begin(spatial_dimension, spatial_dimension); /// Loop on all quadrature points MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & dev_s = *stress_d; Matrix & h = *history_int; /// Compute the first invariant of strain Real Theta = grad_u.trace(); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) { - dev_s(i, j) = 2 * this->mu * (.5 * (grad_u(i,j) + grad_u(j,i)) - 1./3. * Theta *(i == j)); + dev_s(i, j) = + 2 * this->mu * + (.5 * (grad_u(i, j) + grad_u(j, i)) - 1. / 3. * Theta * (i == j)); h(i, j) = 0.; } /// Save the deviator of stress ++stress_d; ++history_int; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialStandardLinearSolidDeviatoric::computeStress(ElementType el_type, GhostType ghost_type) { +template +void MaterialStandardLinearSolidDeviatoric::computeStress( + ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real tau = 0.; // if(std::abs(Ev) > std::numeric_limits::epsilon()) tau = eta / Ev; - Array & stress_dev_vect = stress_dev(el_type, ghost_type); + Array & stress_dev_vect = stress_dev(el_type, ghost_type); Array & history_int_vect = history_integral(el_type, ghost_type); - Array::matrix_iterator stress_d = stress_dev_vect.begin(spatial_dimension, spatial_dimension); - Array::matrix_iterator history_int = history_int_vect.begin(spatial_dimension, spatial_dimension); + Array::matrix_iterator stress_d = + stress_dev_vect.begin(spatial_dimension, spatial_dimension); + Array::matrix_iterator history_int = + history_int_vect.begin(spatial_dimension, spatial_dimension); Matrix s(spatial_dimension, spatial_dimension); Real dt = this->model.getTimeStep(); - Real exp_dt_tau = exp( -dt/tau ); - Real exp_dt_tau_2 = exp( -.5*dt/tau ); + Real exp_dt_tau = exp(-dt / tau); + Real exp_dt_tau_2 = exp(-.5 * dt / tau); Matrix epsilon_d(spatial_dimension, spatial_dimension); Matrix epsilon_v(spatial_dimension, spatial_dimension); /// Loop on all quadrature points MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & dev_s = *stress_d; Matrix & h = *history_int; s.clear(); sigma.clear(); /// Compute the first invariant of strain Real gamma_inf = E_inf / this->E; - Real gamma_v = Ev / this->E; + Real gamma_v = Ev / this->E; this->template gradUToEpsilon(grad_u, epsilon_d); Real Theta = epsilon_d.trace(); epsilon_v.eye(Theta / Real(3.)); epsilon_d -= epsilon_v; Matrix U_rond_prim(spatial_dimension, spatial_dimension); U_rond_prim.eye(gamma_inf * this->kpa * Theta); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) { s(i, j) = 2 * this->mu * epsilon_d(i, j); - h(i, j) = exp_dt_tau * h(i, j) + exp_dt_tau_2 * (s(i, j) - dev_s(i, j)); + h(i, j) = exp_dt_tau * h(i, j) + exp_dt_tau_2 * (s(i, j) - dev_s(i, j)); dev_s(i, j) = s(i, j); - sigma(i, j) = U_rond_prim(i,j) + gamma_inf * s(i, j) + gamma_v * h(i, j); + sigma(i, j) = U_rond_prim(i, j) + gamma_inf * s(i, j) + gamma_v * h(i, j); } /// Save the deviator of stress ++stress_d; ++history_int; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; this->updateDissipatedEnergy(el_type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialStandardLinearSolidDeviatoric::updateDissipatedEnergy(ElementType el_type, - GhostType ghost_type) { +template +void MaterialStandardLinearSolidDeviatoric< + spatial_dimension>::updateDissipatedEnergy(ElementType el_type, + GhostType ghost_type) { AKANTU_DEBUG_IN(); // if(ghost_type == _ghost) return 0.; - + Real tau = 0.; tau = eta / Ev; Real * dis_energy = dissipated_energy(el_type, ghost_type).storage(); - Array & stress_dev_vect = stress_dev(el_type, ghost_type); + Array & stress_dev_vect = stress_dev(el_type, ghost_type); Array & history_int_vect = history_integral(el_type, ghost_type); - Array::matrix_iterator stress_d = stress_dev_vect.begin(spatial_dimension, spatial_dimension); - Array::matrix_iterator history_int = history_int_vect.begin(spatial_dimension, spatial_dimension); + Array::matrix_iterator stress_d = + stress_dev_vect.begin(spatial_dimension, spatial_dimension); + Array::matrix_iterator history_int = + history_int_vect.begin(spatial_dimension, spatial_dimension); Matrix q(spatial_dimension, spatial_dimension); Matrix q_rate(spatial_dimension, spatial_dimension); Matrix epsilon_d(spatial_dimension, spatial_dimension); Matrix epsilon_v(spatial_dimension, spatial_dimension); Real dt = this->model.getTimeStep(); - Real gamma_v = Ev / this->E; + Real gamma_v = Ev / this->E; Real alpha = 1. / (2. * this->mu * gamma_v); /// Loop on all quadrature points MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & dev_s = *stress_d; Matrix & h = *history_int; /// Compute the first invariant of strain this->template gradUToEpsilon(grad_u, epsilon_d); Real Theta = epsilon_d.trace(); epsilon_v.eye(Theta / Real(3.)); epsilon_d -= epsilon_v; q.copy(dev_s); q -= h; q *= gamma_v; q_rate.copy(dev_s); q_rate *= gamma_v; q_rate -= q; q_rate /= tau; for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) - *dis_energy += (epsilon_d(i, j) - alpha * q(i, j)) * q_rate(i, j) * dt; - + *dis_energy += (epsilon_d(i, j) - alpha * q(i, j)) * q_rate(i, j) * dt; /// Save the deviator of stress ++stress_d; ++history_int; ++dis_energy; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -Real MaterialStandardLinearSolidDeviatoric::getDissipatedEnergy() const { +template +Real MaterialStandardLinearSolidDeviatoric< + spatial_dimension>::getDissipatedEnergy() const { AKANTU_DEBUG_IN(); Real de = 0.; /// integrate the dissipated energy for each type of elements - for(auto & type : this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { - de += this->fem.integrate(dissipated_energy(type, _not_ghost), type, - _not_ghost, this->element_filter(type, _not_ghost)); + for (auto & type : + this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { + de += + this->fem.integrate(dissipated_energy(type, _not_ghost), type, + _not_ghost, this->element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return de; } /* -------------------------------------------------------------------------- */ -template -Real MaterialStandardLinearSolidDeviatoric::getDissipatedEnergy(ElementType type, UInt index) const { +template +Real MaterialStandardLinearSolidDeviatoric< + spatial_dimension>::getDissipatedEnergy(ElementType type, + UInt index) const { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = this->fem.getNbIntegrationPoints(type); - Array::const_vector_iterator it = this->dissipated_energy(type, _not_ghost).begin(nb_quadrature_points); + Array::const_vector_iterator it = + this->dissipated_energy(type, _not_ghost).begin(nb_quadrature_points); UInt gindex = (this->element_filter(type, _not_ghost))(index); AKANTU_DEBUG_OUT(); return this->fem.integrate(it[index], type, gindex); } /* -------------------------------------------------------------------------- */ -template -Real MaterialStandardLinearSolidDeviatoric::getEnergy(std::string type) { - if(type == "dissipated") return getDissipatedEnergy(); - else if(type == "dissipated_sls_deviatoric") return getDissipatedEnergy(); - else return MaterialElastic::getEnergy(type); +template +Real MaterialStandardLinearSolidDeviatoric::getEnergy( + const std::string & type) { + if (type == "dissipated") + return getDissipatedEnergy(); + else if (type == "dissipated_sls_deviatoric") + return getDissipatedEnergy(); + else + return MaterialElastic::getEnergy(type); } /* -------------------------------------------------------------------------- */ -template -Real MaterialStandardLinearSolidDeviatoric::getEnergy(std::string energy_id, ElementType type, UInt index) { - if(energy_id == "dissipated") return getDissipatedEnergy(type, index); - else if(energy_id == "dissipated_sls_deviatoric") return getDissipatedEnergy(type, index); - else return MaterialElastic::getEnergy(energy_id, type, index); +template +Real MaterialStandardLinearSolidDeviatoric::getEnergy( + const std::string & energy_id, ElementType type, UInt index) { + if (energy_id == "dissipated") + return getDissipatedEnergy(type, index); + else if (energy_id == "dissipated_sls_deviatoric") + return getDissipatedEnergy(type, index); + else + return MaterialElastic::getEnergy(energy_id, type, + index); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(sls_deviatoric, MaterialStandardLinearSolidDeviatoric); -} // akantu +} // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh b/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh index fa82de98a..a5a70a7ff 100644 --- a/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh +++ b/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh @@ -1,135 +1,136 @@ /** * @file material_standard_linear_solid_deviatoric.hh * * @author David Simon Kammer * @author Nicolas Richart * @author Vladislav Yastrebov * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Oct 08 2015 * * @brief Material Visco-elastic, based on Standard Solid rheological model, * see * [] J.C. Simo, T.J.R. Hughes, "Computational Inelasticity", Springer (1998), * see Sections 10.2 and 10.3 * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_STANDARD_LINEAR_SOLID_DEVIATORIC_HH__ #define __AKANTU_MATERIAL_STANDARD_LINEAR_SOLID_DEVIATORIC_HH__ namespace akantu { /** * Material standard linear solid deviatoric * * * @verbatim E_\inf ------|\/\/\|------ | | ---| |--- | | ----|\/\/\|--[|---- E_v \eta @endverbatim * * keyword : sls_deviatoric * * parameters in the material files : * - E : Initial Young's modulus @f$ E = E_i + E_v @f$ * - eta : viscosity * - Ev : stiffness of the viscous element */ template class MaterialStandardLinearSolidDeviatoric : public MaterialElastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialStandardLinearSolidDeviatoric(SolidMechanicsModel & model, const ID & id = ""); ~MaterialStandardLinearSolidDeviatoric() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter void initMaterial() override; /// update the internal parameters (for modifiable parameters) void updateInternalParameters() override; /// set material to steady state void setToSteadyState(ElementType el_type, GhostType ghost_type = _not_ghost) override; /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; protected: /// update the dissipated energy, is called after the stress have been /// computed void updateDissipatedEnergy(ElementType el_type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// give the dissipated energy for the time step Real getDissipatedEnergy() const; Real getDissipatedEnergy(ElementType type, UInt index) const; /// get the energy using an energy type string for the time step - Real getEnergy(std::string type) override; - Real getEnergy(std::string energy_id, ElementType type, UInt index) override; + Real getEnergy(const std::string & type) override; + Real getEnergy(const std::string & energy_id, ElementType type, + UInt index) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// viscosity, viscous elastic modulus Real eta, Ev, E_inf; /// history of deviatoric stress InternalField stress_dev; /// Internal variable: history integral InternalField history_integral; /// Dissipated energy InternalField dissipated_energy; }; } // namespace akantu #endif /* __AKANTU_MATERIAL_STANDARD_LINEAR_SOLID_DEVIATORIC_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc index 95f999360..de473fdf3 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.cc @@ -1,582 +1,581 @@ /** * @file material_cohesive.cc * * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Jan 12 2016 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "aka_random_generator.hh" #include "dof_synchronizer.hh" #include "fe_engine_template.hh" #include "integrator_gauss.hh" #include "shape_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model, id), facet_filter("facet_filter", id, this->getMemoryID()), fem_cohesive( model.getFEEngineClass("CohesiveFEEngine")), reversible_energy("reversible_energy", *this), total_energy("total_energy", *this), opening("opening", *this), opening_old("opening (old)", *this), tractions("tractions", *this), tractions_old("tractions (old)", *this), contact_tractions("contact_tractions", *this), contact_opening("contact_opening", *this), delta_max("delta max", *this), use_previous_delta_max(false), use_previous_opening(false), damage("damage", *this), sigma_c("sigma_c", *this), normal(0, spatial_dimension, "normal") { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable, "Critical stress"); this->registerParam("delta_c", delta_c, Real(0.), _pat_parsable | _pat_readable, "Critical displacement"); this->element_filter.initialize(this->model->getMesh(), _spatial_dimension = spatial_dimension, _element_kind = _ek_cohesive); // this->model->getMesh().initElementTypeMapArray( // this->element_filter, 1, spatial_dimension, false, _ek_cohesive); if (this->model->getIsExtrinsic()) this->facet_filter.initialize(this->model->getMeshFacets(), _spatial_dimension = spatial_dimension - 1, _element_kind = _ek_regular); // this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1, // spatial_dimension - // 1); this->reversible_energy.initialize(1); this->total_energy.initialize(1); this->tractions_old.initialize(spatial_dimension); this->tractions.initialize(spatial_dimension); this->opening_old.initialize(spatial_dimension); this->contact_tractions.initialize(spatial_dimension); this->contact_opening.initialize(spatial_dimension); this->opening.initialize(spatial_dimension); this->delta_max.initialize(1); this->damage.initialize(1); if (this->model->getIsExtrinsic()) this->sigma_c.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); if (this->use_previous_delta_max) this->delta_max.initializeHistory(); if (this->use_previous_opening) this->opening.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Tractions", tractions); #endif auto & internal_force = const_cast &>(model->getInternalForce()); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) continue; const auto & shapes = fem_cohesive.getShapes(type, ghost_type); auto & traction = tractions(type, ghost_type); auto & contact_traction = contact_tractions(type, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); /// compute @f$t_i N_a@f$ Array * traction_cpy = new Array( nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); auto traction_it = traction.begin(spatial_dimension, 1); auto contact_traction_it = contact_traction.begin(spatial_dimension, 1); auto shapes_filtered_begin = shapes.begin(1, size_of_shapes); auto traction_cpy_it = traction_cpy->begin(spatial_dimension, size_of_shapes); Matrix traction_tmp(spatial_dimension, 1); for (UInt el = 0; el < nb_element; ++el) { UInt current_quad = elem_filter(el) * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++contact_traction_it, ++current_quad, ++traction_cpy_it) { const Matrix & shapes_filtered = shapes_filtered_begin[current_quad]; traction_tmp.copy(*traction_it); traction_tmp += *contact_traction_it; traction_cpy_it->mul(traction_tmp, shapes_filtered); } } /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ Array * int_t_N = new Array( nb_element, spatial_dimension * size_of_shapes, "int_t_N"); fem_cohesive.integrate(*traction_cpy, *int_t_N, spatial_dimension * size_of_shapes, type, ghost_type, elem_filter); delete traction_cpy; int_t_N->extendComponentsInterlaced(2, int_t_N->getNbComponent()); Real * int_t_N_val = int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n) int_t_N_val[n] *= -1.; int_t_N_val += nb_nodes_per_element * spatial_dimension; } /// assemble model->getDOFManager().assembleElementalArrayLocalArray( *int_t_N, internal_force, type, ghost_type, -1, elem_filter); delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { UInt nb_quadrature_points = fem_cohesive.getNbIntegrationPoints(type, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & shapes = fem_cohesive.getShapes(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (!nb_element) continue; UInt size_of_shapes = shapes.getNbComponent(); Array * shapes_filtered = new Array( nb_element * nb_quadrature_points, size_of_shapes, "filtered shapes"); - Real * shapes_val = shapes.storage(); Real * shapes_filtered_val = shapes_filtered->storage(); UInt * elem_filter_val = elem_filter.storage(); for (UInt el = 0; el < nb_element; ++el) { - shapes_val = shapes.storage() + + auto shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } /** * compute A matrix @f$ \mathbf{A} = \left[\begin{array}{c c c c c c c c c c *c c} * 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ * 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\ * 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ * 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\ * 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\ * 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 * \end{array} \right]@f$ **/ // UInt size_of_A = // spatial_dimension*size_of_shapes*spatial_dimension*nb_nodes_per_element; // Real * A = new Real[size_of_A]; // memset(A, 0, size_of_A*sizeof(Real)); Matrix A(spatial_dimension * size_of_shapes, spatial_dimension * nb_nodes_per_element); for (UInt i = 0; i < spatial_dimension * size_of_shapes; ++i) { A(i, i) = 1; A(i, i + spatial_dimension * size_of_shapes) = -1; } // compute traction. This call is not necessary for the linear // cohesive law that, currently, is the only one used for the // extrinsic approach. // if (!model->getIsExtrinsic()){ // computeTraction(ghost_type); //} /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} /// @f$ Array * tangent_stiffness_matrix = new Array( nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix"); // Array * normal = new Array(nb_element * // nb_quadrature_points, spatial_dimension, "normal"); normal.resize(nb_quadrature_points); computeNormal(model->getCurrentPosition(), normal, type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ // computeOpening(model->getDisplacement(), opening(type, ghost_type), type, // ghost_type); tangent_stiffness_matrix->clear(); computeTangentTraction(type, *tangent_stiffness_matrix, normal, ghost_type); // delete normal; UInt size_at_nt_d_n_a = spatial_dimension * nb_nodes_per_element * spatial_dimension * nb_nodes_per_element; Array * at_nt_d_n_a = new Array( nb_element * nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Array::iterator> shapes_filt_it = shapes_filtered->begin(size_of_shapes); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Array::matrix_iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Array::matrix_iterator At_Nt_D_N_A_end = at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Matrix N(spatial_dimension, spatial_dimension * size_of_shapes); Matrix N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); Matrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for (; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { N.clear(); /** * store the shapes in voigt notations matrix @f$\mathbf{N} = * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$ **/ for (UInt i = 0; i < spatial_dimension; ++i) for (UInt n = 0; n < size_of_shapes; ++n) N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n); /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} *{\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ N_A.mul(N, A); D_N_A.mul(*D_it, N_A); (*At_Nt_D_N_A_it).mul(D_N_A, N_A); } delete tangent_stiffness_matrix; delete shapes_filtered; Array * K_e = new Array(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive.integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, type, ghost_type, elem_filter); delete at_nt_d_n_a; model->getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *K_e, type, ghost_type, _unsymmetric, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Openings", opening); #endif for (auto & type : element_filter.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { Array & elem_filter = element_filter(type, ghost_type); UInt nb_element = elem_filter.size(); if (nb_element == 0) continue; UInt nb_quadrature_points = nb_element * fem_cohesive.getNbIntegrationPoints(type, ghost_type); normal.resize(nb_quadrature_points); /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normal, type, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(type, ghost_type), type, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normal, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); if (type == _cohesive_1d_2) fem_cohesive.computeNormalsOnIntegrationPoints(position, normal, type, ghost_type); else { #define COMPUTE_NORMAL(type) \ fem_cohesive.getShapeFunctions() \ .computeNormalsOnIntegrationPoints( \ position, normal, ghost_type, element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem_cohesive = this->model->getFEEngineClass("CohesiveFEEngine"); #define COMPUTE_OPENING(type) \ fem_cohesive.getShapeFunctions() \ .interpolateOnIntegrationPoints( \ displacement, opening, spatial_dimension, ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeEnergies() { AKANTU_DEBUG_IN(); Vector b(spatial_dimension); Vector h(spatial_dimension); for (auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost, _ek_cohesive)) { auto erev = reversible_energy(type, _not_ghost).begin(); auto etot = total_energy(type, _not_ghost).begin(); auto traction_it = tractions(type, _not_ghost).begin(spatial_dimension); auto traction_old_it = tractions_old(type, _not_ghost).begin(spatial_dimension); auto opening_it = opening(type, _not_ghost).begin(spatial_dimension); auto opening_old_it = opening_old(type, _not_ghost).begin(spatial_dimension); auto traction_end = tractions(type, _not_ghost).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration b = *opening_it; b -= *opening_old_it; h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } } /// update old values GhostType ghost_type = _not_ghost; for (auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost, _ek_cohesive)) { tractions_old(type, ghost_type).copy(tractions(type, ghost_type)); opening_old(type, ghost_type).copy(opening(type, ghost_type)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate reversible energy for each type of elements for (auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost, _ek_cohesive)) { erev += fem_cohesive.integrate(reversible_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate dissipated energy for each type of elements for (auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost, _ek_cohesive)) { Array dissipated_energy(total_energy(type, _not_ghost)); dissipated_energy -= reversible_energy(type, _not_ghost); edis += fem_cohesive.integrate(dissipated_energy, type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getContactEnergy() { AKANTU_DEBUG_IN(); Real econ = 0.; /// integrate contact energy for each type of elements for (auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost, _ek_cohesive)) { auto & el_filter = element_filter(type, _not_ghost); UInt nb_quad_per_el = fem_cohesive.getNbIntegrationPoints(type, _not_ghost); UInt nb_quad_points = el_filter.size() * nb_quad_per_el; Array contact_energy(nb_quad_points); auto contact_traction_it = contact_tractions(type, _not_ghost).begin(spatial_dimension); auto contact_opening_it = contact_opening(type, _not_ghost).begin(spatial_dimension); /// loop on each quadrature point for (UInt q = 0; q < nb_quad_points; ++contact_traction_it, ++contact_opening_it, ++q) { contact_energy(q) = .5 * contact_traction_it->dot(*contact_opening_it); } econ += fem_cohesive.integrate(contact_energy, type, _not_ghost, el_filter); } AKANTU_DEBUG_OUT(); return econ; } /* -------------------------------------------------------------------------- */ -Real MaterialCohesive::getEnergy(std::string type) { +Real MaterialCohesive::getEnergy(const std::string & type) { AKANTU_DEBUG_IN(); if (type == "reversible") return getReversibleEnergy(); else if (type == "dissipated") return getDissipatedEnergy(); else if (type == "cohesive contact") return getContactEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh index f33d2e307..4643f98ac 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh @@ -1,256 +1,256 @@ /** * @file material_cohesive.hh * * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Jan 12 2016 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" /* -------------------------------------------------------------------------- */ #include "cohesive_internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_HH__ #define __AKANTU_MATERIAL_COHESIVE_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } namespace akantu { class MaterialCohesive : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineCohesiveType = FEEngineTemplate; public: MaterialCohesive(SolidMechanicsModel & model, const ID & id = ""); ~MaterialCohesive() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter void initMaterial() override; /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// assemble residual void assembleInternalForces(GhostType ghost_type = _not_ghost) override; /// compute reversible and total energies by element void computeEnergies(); /// check stress for cohesive elements' insertion, by default it /// also updates the cohesive elements' data virtual void checkInsertion(bool /*check_only*/ = false) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// check delta_max for cohesive elements in case of no convergence /// in the solveStep (only for extrinsic-implicit) virtual void checkDeltaMax(GhostType /*ghost_type*/ = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// reset variables when convergence is reached (only for /// extrinsic-implicit) virtual void resetVariables(GhostType /*ghost_type*/ = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// interpolate stress on given positions for each element (empty /// implemantation to avoid the generic call to be done on cohesive elements) virtual void interpolateStress(const ElementType /*type*/, Array & /*result*/){}; /// compute the stresses void computeAllStresses(GhostType /*ghost_type*/ = _not_ghost) override{}; // add the facet to be handled by the material UInt addFacet(const Element & element); protected: virtual void computeTangentTraction(const ElementType & /*el_type*/, Array & /*tangent_matrix*/, const Array & /*normal*/, GhostType /*ghost_type*/ = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the normal void computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type); /// compute the opening void computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type); template void computeNormal(const Array & position, Array & normal, GhostType ghost_type); /// assemble stiffness void assembleStiffnessMatrix(GhostType ghost_type) override; /// constitutive law virtual void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// parallelism functions inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening, Real); /// get the traction AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Traction, tractions, Real); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /// get facet filter AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO(FacetFilter, facet_filter, const ElementTypeMapArray &); // AKANTU_GET_MACRO(ElementFilter, element_filter, const // ElementTypeMapArray &); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// compute contact energy Real getContactEnergy(); /// get energy - Real getEnergy(std::string type) override; + Real getEnergy(const std::string &type) override; /// return the energy (identified by id) for the provided element - Real getEnergy(std::string energy_id, ElementType type, UInt index) override { + Real getEnergy(const std::string &energy_id, ElementType type, UInt index) override { return Material::getEnergy(energy_id, type, index); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// list of facets assigned to this material ElementTypeMapArray facet_filter; /// Link to the cohesive fem object in the model FEEngine & fem_cohesive; private: /// reversible energy by quadrature point CohesiveInternalField reversible_energy; /// total energy by quadrature point CohesiveInternalField total_energy; protected: /// opening in all elements and quadrature points CohesiveInternalField opening; /// opening in all elements and quadrature points (previous time step) CohesiveInternalField opening_old; /// traction in all elements and quadrature points CohesiveInternalField tractions; /// traction in all elements and quadrature points (previous time step) CohesiveInternalField tractions_old; /// traction due to contact CohesiveInternalField contact_tractions; /// normal openings for contact tractions CohesiveInternalField contact_opening; /// maximum displacement CohesiveInternalField delta_max; /// tell if the previous delta_max state is needed (in iterative schemes) bool use_previous_delta_max; /// tell if the previous opening state is needed (in iterative schemes) bool use_previous_opening; /// damage CohesiveInternalField damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; /// critical stress RandomInternalField sigma_c; /// critical displacement Real delta_c; /// array to temporarily store the normals Array normal; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_inline_impl.cc" } // namespace akantu #include "cohesive_internal_field_tmpl.hh" #endif /* __AKANTU_MATERIAL_COHESIVE_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc index 78d1cfaa1..90e386726 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_parallel.cc @@ -1,397 +1,400 @@ /** * @file solid_mechanics_model_cohesive_parallel.cc * * @author Marco Vocialta * * * @brief Functions for parallel cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "element_synchronizer.hh" #include "solid_mechanics_model_cohesive.hh" #include "solid_mechanics_model_tmpl.hh" #include "communicator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::synchronizeGhostFacetsConnectivity() { AKANTU_DEBUG_IN(); const Communicator & comm = mesh.getCommunicator(); Int psize = comm.getNbProc(); if (psize > 1) { /// get global connectivity for not ghost facets global_connectivity = new ElementTypeMapArray("global_connectivity", id); auto & mesh_facets = inserter->getMeshFacets(); global_connectivity->initialize( mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true, _with_nb_nodes_per_element = true, _element_kind = _ek_regular, _ghost_type = _not_ghost); mesh_facets.getGlobalConnectivity(*global_connectivity); /// communicate synchronize(_gst_smmc_facets_conn); /// flip facets MeshUtils::flipFacets(mesh_facets, *global_connectivity, _ghost); delete global_connectivity; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateCohesiveSynchronizers() { /// update synchronizers if needed if (not mesh.isDistributed()) return; auto & mesh_facets = inserter->getMeshFacets(); auto & facet_synchronizer = mesh_facets.getElementSynchronizer(); const auto & cfacet_synchronizer = facet_synchronizer; // update the cohesive element synchronizer cohesive_synchronizer->updateSchemes( [&](auto && scheme, auto && proc, auto && direction) { auto & facet_scheme = cfacet_synchronizer.getCommunications().getScheme(proc, direction); for (auto && facet : facet_scheme) { const auto & connected_elements = mesh_facets.getElementToSubelement( facet.type, facet.ghost_type)(facet.element); // slow access here const auto & cohesive_element = connected_elements[1]; auto && cohesive_type = FEEngine::getCohesiveElementType(facet.type); auto old_nb_cohesive_elements = mesh.getNbElement(cohesive_type, facet.ghost_type); old_nb_cohesive_elements -= mesh.getData("facet_to_double", facet.type, facet.ghost_type).size(); if (cohesive_element.kind() == _ek_cohesive and cohesive_element.element >= old_nb_cohesive_elements) { scheme.push_back(cohesive_element); } } }); const auto & comm = mesh.getCommunicator(); auto && my_rank = comm.whoAmI(); // update the facet stress synchronizer facet_stress_synchronizer->updateSchemes([&](auto && scheme, auto && proc, auto && /*direction*/) { auto it_element = scheme.begin(); for (auto && element : scheme) { auto && facet_check = inserter->getCheckFacets( element.type, element.ghost_type)(element.element); // slow access // here if (facet_check) { auto && connected_elements = mesh_facets.getElementToSubelement( element.type, element.ghost_type)(element.element); // slow access // here auto && rank_left = facet_synchronizer.getRank(connected_elements[0]); auto && rank_right = facet_synchronizer.getRank(connected_elements[1]); // keep element if the element is still a boundary element between two // processors if ((rank_left == Int(proc) and rank_right == my_rank) or (rank_left == my_rank and rank_right == Int(proc))) { *it_element = element; ++it_element; } } } scheme.resize(it_element - scheme.begin()); }); } /* -------------------------------------------------------------------------- */ template void SolidMechanicsModelCohesive::packFacetStressDataHelper( const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements) const { packUnpackFacetStressDataHelper( const_cast &>(data_to_pack), buffer, elements); } /* -------------------------------------------------------------------------- */ template void SolidMechanicsModelCohesive::unpackFacetStressDataHelper( ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements) const { packUnpackFacetStressDataHelper(data_to_unpack, buffer, elements); } /* -------------------------------------------------------------------------- */ template void SolidMechanicsModelCohesive::packUnpackFacetStressDataHelper( ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements) const { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt nb_quad_per_elem = 0; UInt sp2 = spatial_dimension * spatial_dimension; UInt nb_component = sp2 * 2; bool element_rank = false; Mesh & mesh_facets = inserter->getMeshFacets(); Array * vect = nullptr; Array> * element_to_facet = nullptr; auto & fe_engine = this->getFEEngine("FacetsFEEngine"); for (auto && el : elements) { + if (el.type == _not_defined) + AKANTU_EXCEPTION("packUnpackFacetStressDataHelper called with wrong inputs"); + if (el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; vect = &data_to_pack(el.type, el.ghost_type); element_to_facet = &(mesh_facets.getElementToSubelement(el.type, el.ghost_type)); nb_quad_per_elem = fe_engine.getNbIntegrationPoints(el.type, el.ghost_type); } if (pack_helper) element_rank = (*element_to_facet)(el.element)[0].ghost_type != _not_ghost; else element_rank = (*element_to_facet)(el.element)[0].ghost_type == _not_ghost; for (UInt q = 0; q < nb_quad_per_elem; ++q) { Vector data(vect->storage() + (el.element * nb_quad_per_elem + q) * nb_component + element_rank * sp2, sp2); if (pack_helper) buffer << data; else buffer >> data; } } } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::getNbQuadsForFacetCheck( const Array & elements) const { UInt nb_quads = 0; UInt nb_quad_per_facet = 0; ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; auto & fe_engine = this->getFEEngine("FacetsFEEngine"); for(auto & el : elements) { if (el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; nb_quad_per_facet = fe_engine .getNbIntegrationPoints(el.type, el.ghost_type); } nb_quads += nb_quad_per_facet; } return nb_quads; } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::getNbData( const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; if (elements.size() == 0) return 0; /// regular element case if (elements(0).kind() == _ek_regular) { switch (tag) { case _gst_smmc_facets: { size += elements.size() * sizeof(bool); break; } case _gst_smmc_facets_conn: { UInt nb_nodes = Mesh::getNbNodesPerElementList(elements); size += nb_nodes * sizeof(UInt); break; } case _gst_smmc_facets_stress: { UInt nb_quads = getNbQuadsForFacetCheck(elements); size += nb_quads * spatial_dimension * spatial_dimension * sizeof(Real); break; } default: { size += SolidMechanicsModel::getNbData(elements, tag); } } } /// cohesive element case else if (elements(0).kind() == _ek_cohesive) { switch (tag) { case _gst_material_id: { size += elements.size() * sizeof(UInt); break; } case _gst_smm_boundary: { UInt nb_nodes_per_element = 0; for (auto && el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } // force, displacement, boundary size += nb_nodes_per_element * spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } default: break; } if (tag != _gst_material_id && tag != _gst_smmc_facets) { splitByMaterial(elements, [&](auto && mat, auto && elements) { size += mat.getNbData(elements, tag); }); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::packData( CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); if (elements.size() == 0) return; if (elements(0).kind() == _ek_regular) { switch (tag) { case _gst_smmc_facets: { packElementalDataHelper(inserter->getInsertionFacetsByElement(), buffer, elements, false, getFEEngine()); break; } case _gst_smmc_facets_conn: { packElementalDataHelper(*global_connectivity, buffer, elements, false, getFEEngine()); break; } case _gst_smmc_facets_stress: { packFacetStressDataHelper(facet_stress, buffer, elements); break; } default: { SolidMechanicsModel::packData(buffer, elements, tag); } } } else if (elements(0).kind() == _ek_cohesive) { switch (tag) { case _gst_material_id: { packElementalDataHelper(material_index, buffer, elements, false, getFEEngine("CohesiveFEEngine")); break; } case _gst_smm_boundary: { packNodalDataHelper(*internal_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != _gst_material_id && tag != _gst_smmc_facets) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.packData(buffer, elements, tag); }); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); if (elements.size() == 0) return; if (elements(0).kind() == _ek_regular) { switch (tag) { case _gst_smmc_facets: { unpackElementalDataHelper(inserter->getInsertionFacetsByElement(), buffer, elements, false, getFEEngine()); break; } case _gst_smmc_facets_conn: { unpackElementalDataHelper(*global_connectivity, buffer, elements, false, getFEEngine()); break; } case _gst_smmc_facets_stress: { unpackFacetStressDataHelper(facet_stress, buffer, elements); break; } default: { SolidMechanicsModel::unpackData(buffer, elements, tag); } } } else if (elements(0).kind() == _ek_cohesive) { switch (tag) { case _gst_material_id: { unpackElementalDataHelper(material_index, buffer, elements, false, getFEEngine("CohesiveFEEngine")); break; } case _gst_smm_boundary: { unpackNodalDataHelper(*internal_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != _gst_material_id && tag != _gst_smmc_facets) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.unpackData(buffer, elements, tag); }); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/solver/sparse_solver.cc b/src/solver/sparse_solver.cc index e0060e98b..47a4dbb09 100644 --- a/src/solver/sparse_solver.cc +++ b/src/solver/sparse_solver.cc @@ -1,86 +1,86 @@ /** * @file solver.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Tue Jan 19 2016 * * @brief Solver interface class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh.hh" #include "dof_manager.hh" #include "sparse_solver.hh" #include "communicator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SparseSolver::SparseSolver(DOFManager & dof_manager, const ID & matrix_id, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), Parsable(_st_solver, id), _dof_manager(dof_manager), matrix_id(matrix_id), communicator(dof_manager.getCommunicator()) { AKANTU_DEBUG_IN(); // OK this is fishy... this->communicator.registerEventHandler(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseSolver::~SparseSolver() { AKANTU_DEBUG_IN(); - this->destroyInternalData(); + //this->destroyInternalData(); this->communicator.unregisterEventHandler(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolver::beforeStaticSolverDestroy() { AKANTU_DEBUG_IN(); try { this->destroyInternalData(); } catch (...) { } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolver::createSynchronizerRegistry() { // this->synch_registry = new SynchronizerRegistry(this); } void SparseSolver::onCommunicatorFinalize() { this->destroyInternalData(); } } // akantu diff --git a/src/solver/sparse_solver_mumps.cc b/src/solver/sparse_solver_mumps.cc index b84b588af..f2a1e9a7e 100644 --- a/src/solver/sparse_solver_mumps.cc +++ b/src/solver/sparse_solver_mumps.cc @@ -1,440 +1,447 @@ /** * @file sparse_solver_mumps.cc * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Tue Jan 19 2016 * * @brief implem of SparseSolverMumps class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @subsection Ctrl_param Control parameters * * ICNTL(1), * ICNTL(2), * ICNTL(3) : output streams for error, diagnostics, and global messages * * ICNTL(4) : verbose level : 0 no message - 4 all messages * * ICNTL(5) : type of matrix, 0 assembled, 1 elementary * * ICNTL(6) : control the permutation and scaling(default 7) see mumps doc for * more information * * ICNTL(7) : determine the pivot order (default 7) see mumps doc for more * information * * ICNTL(8) : describe the scaling method used * * ICNTL(9) : 1 solve A x = b, 0 solve At x = b * * ICNTL(10) : number of iterative refinement when NRHS = 1 * * ICNTL(11) : > 0 return statistics * * ICNTL(12) : only used for SYM = 2, ordering strategy * * ICNTL(13) : * * ICNTL(14) : percentage of increase of the estimated working space * * ICNTL(15-17) : not used * * ICNTL(18) : only used if ICNTL(5) = 0, 0 matrix centralized, 1 structure on * host and mumps give the mapping, 2 structure on host and distributed matrix * for facto, 3 distributed matrix * * ICNTL(19) : > 0, Shur complement returned * * ICNTL(20) : 0 rhs dense, 1 rhs sparse * * ICNTL(21) : 0 solution in rhs, 1 solution distributed in ISOL_loc and SOL_loc * allocated by user * * ICNTL(22) : 0 in-core, 1 out-of-core * * ICNTL(23) : maximum memory allocatable by mumps pre proc * * ICNTL(24) : controls the detection of "null pivot rows" * * ICNTL(25) : * * ICNTL(26) : * * ICNTL(27) : * * ICNTL(28) : 0 automatic choice, 1 sequential analysis, 2 parallel analysis * * ICNTL(29) : 0 automatic choice, 1 PT-Scotch, 2 ParMetis */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "dof_manager_default.hh" #include "dof_synchronizer.hh" #include "sparse_matrix_aij.hh" #if defined(AKANTU_USE_MPI) #include "mpi_communicator_data.hh" #endif #include "sparse_solver_mumps.hh" /* -------------------------------------------------------------------------- */ // static std::ostream & operator <<(std::ostream & stream, const DMUMPS_STRUC_C // & _this) { // stream << "DMUMPS Data [" << std::endl; // stream << " + job : " << _this.job << std::endl; // stream << " + par : " << _this.par << std::endl; // stream << " + sym : " << _this.sym << std::endl; // stream << " + comm_fortran : " << _this.comm_fortran << std::endl; // stream << " + nz : " << _this.nz << std::endl; // stream << " + irn : " << _this.irn << std::endl; // stream << " + jcn : " << _this.jcn << std::endl; // stream << " + nz_loc : " << _this.nz_loc << std::endl; // stream << " + irn_loc : " << _this.irn_loc << std::endl; // stream << " + jcn_loc : " << _this.jcn_loc << std::endl; // stream << "]"; // return stream; // } namespace akantu { /* -------------------------------------------------------------------------- */ SparseSolverMumps::SparseSolverMumps(DOFManagerDefault & dof_manager, const ID & matrix_id, const ID & id, const MemoryID & memory_id) : SparseSolver(dof_manager, matrix_id, id, memory_id), dof_manager(dof_manager), master_rhs_solution(0, 1) { AKANTU_DEBUG_IN(); this->prank = communicator.whoAmI(); #ifdef AKANTU_USE_MPI this->parallel_method = _fully_distributed; #else // AKANTU_USE_MPI this->parallel_method = _not_parallel; #endif // AKANTU_USE_MPI AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseSolverMumps::~SparseSolverMumps() { AKANTU_DEBUG_IN(); + mumpsDataDestroy(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void SparseSolverMumps::destroyInternalData() { +void SparseSolverMumps::mumpsDataDestroy() { if (this->is_initialized) { this->mumps_data.job = _smj_destroy; // destroy dmumps_c(&this->mumps_data); this->is_initialized = false; } } +/* -------------------------------------------------------------------------- */ +void SparseSolverMumps::destroyInternalData() { + mumpsDataDestroy(); +} + /* -------------------------------------------------------------------------- */ void SparseSolverMumps::checkInitialized() { if (this->is_initialized) return; this->initialize(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::setOutputLevel() { // Output setup icntl(1) = 0; // error output icntl(2) = 0; // diagnostics output icntl(3) = 0; // information icntl(4) = 0; #if !defined(AKANTU_NDEBUG) DebugLevel dbg_lvl = debug::debugger.getDebugLevel(); if (AKANTU_DEBUG_TEST(dblDump)) { strcpy(this->mumps_data.write_problem, "mumps_matrix.mtx"); } // clang-format off icntl(1) = (dbg_lvl >= dblWarning) ? 6 : 0; icntl(3) = (dbg_lvl >= dblInfo) ? 6 : 0; icntl(2) = (dbg_lvl >= dblTrace) ? 6 : 0; icntl(4) = dbg_lvl >= dblDump ? 4 : dbg_lvl >= dblTrace ? 3 : dbg_lvl >= dblInfo ? 2 : dbg_lvl >= dblWarning ? 1 : 0; // clang-format on #endif } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::initMumpsData() { auto & A = dof_manager.getMatrix(matrix_id); // Default Scaling icntl(8) = 77; // Assembled matrix icntl(5) = 0; /// Default centralized dense second member icntl(20) = 0; icntl(21) = 0; // automatic choice for analysis icntl(28) = 0; UInt size = A.size(); if (prank == 0) { this->master_rhs_solution.resize(size); } this->mumps_data.nz_alloc = 0; this->mumps_data.n = size; switch (this->parallel_method) { case _fully_distributed: icntl(18) = 3; // fully distributed this->mumps_data.nz_loc = A.getNbNonZero(); this->mumps_data.irn_loc = A.getIRN().storage(); this->mumps_data.jcn_loc = A.getJCN().storage(); break; case _not_parallel: case _master_slave_distributed: icntl(18) = 0; // centralized if (prank == 0) { this->mumps_data.nz = A.getNbNonZero(); this->mumps_data.irn = A.getIRN().storage(); this->mumps_data.jcn = A.getJCN().storage(); } else { this->mumps_data.nz = 0; this->mumps_data.irn = nullptr; this->mumps_data.jcn = nullptr; } break; default: AKANTU_DEBUG_ERROR("This case should not happen!!"); } } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::initialize() { AKANTU_DEBUG_IN(); this->mumps_data.par = 1; // The host is part of computations switch (this->parallel_method) { case _not_parallel: break; case _master_slave_distributed: this->mumps_data.par = 0; // The host is not part of the computations /* FALLTHRU */ /* [[fallthrough]]; un-comment when compiler will get it */ case _fully_distributed: #ifdef AKANTU_USE_MPI const auto & mpi_data = dynamic_cast( communicator.getCommunicatorData()); MPI_Comm mpi_comm = mpi_data.getMPICommunicator(); this->mumps_data.comm_fortran = MPI_Comm_c2f(mpi_comm); #else AKANTU_DEBUG_ERROR( "You cannot use parallel method to solve without activating MPI"); #endif break; } const auto & A = dof_manager.getMatrix(matrix_id); this->mumps_data.sym = 2 * (A.getMatrixType() == _symmetric); this->prank = communicator.whoAmI(); this->setOutputLevel(); this->mumps_data.job = _smj_initialize; // initialize dmumps_c(&this->mumps_data); this->setOutputLevel(); this->is_initialized = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::analysis() { AKANTU_DEBUG_IN(); initMumpsData(); this->mumps_data.job = _smj_analyze; // analyze dmumps_c(&this->mumps_data); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::factorize() { AKANTU_DEBUG_IN(); auto & A = dof_manager.getMatrix(matrix_id); if (parallel_method == _fully_distributed) this->mumps_data.a_loc = A.getA().storage(); else { if (prank == 0) this->mumps_data.a = A.getA().storage(); } this->mumps_data.job = _smj_factorize; // factorize dmumps_c(&this->mumps_data); this->printError(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::solve(Array & x, const Array & b) { auto & synch = this->dof_manager.getSynchronizer(); if (this->prank == 0) { this->master_rhs_solution.resize(this->dof_manager.getSystemSize()); synch.gather(b, this->master_rhs_solution); } else { synch.gather(b); } this->solveInternal(); if (this->prank == 0) { synch.scatter(x, this->master_rhs_solution); } else { synch.scatter(x); } } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::solve() { this->master_rhs_solution.copy(this->dof_manager.getGlobalResidual()); this->solveInternal(); this->dof_manager.setGlobalSolution(this->master_rhs_solution); this->dof_manager.splitSolutionPerDOFs(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::solveInternal() { AKANTU_DEBUG_IN(); this->checkInitialized(); const auto & A = dof_manager.getMatrix(matrix_id); this->setOutputLevel(); if (this->last_profile_release != A.getProfileRelease()) { this->analysis(); this->last_profile_release = A.getProfileRelease(); } if (AKANTU_DEBUG_TEST(dblDump)) { std::stringstream sstr; sstr << prank << ".mtx"; A.saveMatrix("solver_mumps" + sstr.str()); } if (this->last_value_release != A.getValueRelease()) { this->factorize(); this->last_value_release = A.getValueRelease(); } if (prank == 0) { this->mumps_data.rhs = this->master_rhs_solution.storage(); } this->mumps_data.job = _smj_solve; // solve dmumps_c(&this->mumps_data); this->printError(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::printError() { Vector _info_v(2); _info_v[0] = info(1); // to get errors _info_v[1] = -info(1); // to get warnings dof_manager.getCommunicator().allReduce(_info_v, SynchronizerOperation::_min); _info_v[1] = -_info_v[1]; if (_info_v[0] < 0) { // < 0 is an error switch (_info_v[0]) { case -10: AKANTU_DEBUG_ERROR("The matrix is singular"); break; case -9: { icntl(14) += 10; if (icntl(14) != 90) { // std::cout << "Dynamic memory increase of 10%" << std::endl; AKANTU_DEBUG_WARNING("MUMPS dynamic memory is insufficient it will be " "increased allowed to use 10% more"); // change releases to force a recompute this->last_value_release--; this->last_profile_release--; this->solve(); } else { AKANTU_DEBUG_ERROR("The MUMPS workarray is too small INFO(2)=" << info(2) << "No further increase possible"); } break; } default: AKANTU_DEBUG_ERROR("Error in mumps during solve process, check mumps " "user guide INFO(1) = " << _info_v[1]); } } else if (_info_v[1] > 0) { AKANTU_DEBUG_WARNING("Warning in mumps during solve process, check mumps " "user guide INFO(1) = " << _info_v[1]); } } } // akantu diff --git a/src/solver/sparse_solver_mumps.hh b/src/solver/sparse_solver_mumps.hh index f200333b5..6be46c20c 100644 --- a/src/solver/sparse_solver_mumps.hh +++ b/src/solver/sparse_solver_mumps.hh @@ -1,157 +1,160 @@ /** * @file sparse_solver_mumps.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Jan 19 2016 * * @brief Solver class implementation for the mumps solver * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "sparse_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_MUMPS_HH__ #define __AKANTU_SOLVER_MUMPS_HH__ namespace akantu { class DOFManagerDefault; class SparseMatrixAIJ; } namespace akantu { class SparseSolverMumps : public SparseSolver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SparseSolverMumps(DOFManagerDefault & dof_manager, const ID & matrix_id, const ID & id = "sparse_solver_mumps", const MemoryID & memory_id = 0); ~SparseSolverMumps() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build the profile and do the analysis part void initialize() override; /// analysis (symbolic facto + permutations) void analysis() override; /// factorize the matrix void factorize() override; /// solve the system virtual void solve(Array & x, const Array & b); /// solve using residual and solution from the dof_manager void solve() override; private: /// print the error if any happened in mumps void printError(); /// solve the system with master_rhs_solution as b and x void solveInternal(); /// set internal values; void initMumpsData(); /// set the level of verbosity of mumps based on the debug level of akantu void setOutputLevel(); protected: /// de-initialize the internal data void destroyInternalData() override; /// check if initialized and except if it is not the case void checkInitialized(); +private: + + void mumpsDataDestroy(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ private: /// access the control variable inline Int & icntl(UInt i) { return mumps_data.icntl[i - 1]; } /// access the results info inline Int & info(UInt i) { return mumps_data.info[i - 1]; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// DOFManager used by the Mumps implementation of the SparseSolver DOFManagerDefault & dof_manager; /// Full right hand side on the master processors and solution after solve Array master_rhs_solution; /// mumps data DMUMPS_STRUC_C mumps_data; /// Rank of the current process UInt prank; /// matrix release at last solve UInt last_profile_release{UInt(-1)}; /// matrix release at last solve UInt last_value_release{UInt(-1)}; /// check if the solver data are initialized bool is_initialized{false}; /* ------------------------------------------------------------------------ */ /* Local types */ /* ------------------------------------------------------------------------ */ private: SolverParallelMethod parallel_method; // bool rhs_is_local; enum SolverMumpsJob { _smj_initialize = -1, _smj_analyze = 1, _smj_factorize = 2, _smj_solve = 3, _smj_analyze_factorize = 4, _smj_factorize_solve = 5, _smj_complete = 6, // analyze, factorize, solve _smj_destroy = -2 }; }; } // akantu #endif /* __AKANTU_SOLVER_MUMPS_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc index dc4b498f2..0e4ea4526 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc @@ -1,351 +1,351 @@ /** * @file test_cohesive_intrinsic_tetrahedron.cc * * @author Marco Vocialta * * @date creation: Tue Aug 27 2013 * @date last modification: Thu Oct 15 2015 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; class Checker { public: Checker(const SolidMechanicsModelCohesive & model, const Array & elements, ElementType type); void check(const Vector & opening, const Matrix & rotation) { checkTractions(opening, rotation); checkEquilibrium(); computeEnergy(opening); } void updateDisplacement(const Vector & increment); protected: void checkTractions(const Vector & opening, const Matrix & rotation); void checkEquilibrium(); void checkResidual(const Matrix & rotation); void computeEnergy(const Vector & opening); private: std::set nodes_to_check; const SolidMechanicsModelCohesive & model; ElementType type; - const Array & elements; + //const Array & elements; const Material & mat_cohesive; Real sigma_c; const Real beta; const Real G_c; const Real delta_0; const Real kappa; Real delta_c; const UInt spatial_dimension; const ElementType type_facet; const ElementType type_cohesive; const Array & traction; const Array & damage; const UInt nb_quad_per_el; const UInt nb_element; const Real beta2_kappa2; const Real beta2_kappa; Vector theoretical_traction; Vector traction_old; Vector opening_old; Real Ed; }; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("material_tetrahedron.dat", argc, argv); // debug::setDebugLevel(dblDump); const UInt spatial_dimension = 3; const UInt max_steps = 60; const Real increment_constant = 0.01; Math::setTolerance(1.e-12); const ElementType type = _tetrahedron_10; Mesh mesh(spatial_dimension); mesh.read("tetrahedron.msh"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull(); model.limitInsertion(_x, -0.01, 0.01); model.insertIntrinsicElements(); Array & boundary = model.getBlockedDOFs(); boundary.set(true); UInt nb_element = mesh.getNbElement(type); model.setBaseName("intrinsic_tetrahedron"); model.addDumpFieldVector("displacement"); model.addDumpField("internal_force"); model.dump(); model.setBaseNameToDumper("cohesive elements", "cohesive_elements_tetrahedron"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); model.dump("cohesive elements"); /// find elements to displace Array elements; Vector bary(spatial_dimension); for (UInt el = 0; el < nb_element; ++el) { mesh.getBarycenter(el, type, bary.storage()); if (bary(_x) > 0.01) elements.push_back(el); } /// rotate mesh Real angle = 1.; // clang-format off Matrix rotation{ {std::cos(angle), std::sin(angle) * -1., 0.}, {std::sin(angle), std::cos(angle), 0.}, {0., 0., 1.}}; // clang-format on Vector increment_tmp{increment_constant, 2. * increment_constant, 3. * increment_constant}; Vector increment = rotation * increment_tmp; auto & position = mesh.getNodes(); auto position_it = position.begin(spatial_dimension); auto position_end = position.end(spatial_dimension); for (; position_it != position_end; ++position_it) { auto & pos = *position_it; pos = rotation * pos; } model.dump(); model.dump("cohesive elements"); /// find nodes to check Checker checker(model, elements, type); checker.updateDisplacement(increment); Real theoretical_Ed = 0; Vector opening(spatial_dimension, 0.); Vector opening_old(spatial_dimension, 0.); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { model.solveStep(); model.dump(); model.dump("cohesive elements"); opening += increment_tmp; checker.check(opening, rotation); checker.updateDisplacement(increment); } model.dump(); model.dump("cohesive elements"); Real Ed = model.getEnergy("dissipated"); theoretical_Ed *= 4.; std::cout << Ed << " " << theoretical_Ed << std::endl; if (!Math::are_float_equal(Ed, theoretical_Ed) || std::isnan(Ed)) { std::cout << "The dissipated energy is incorrect" << std::endl; finalize(); return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_intrinsic_tetrahedron was passed!" << std::endl; return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void Checker::updateDisplacement(const Vector & increment) { Mesh & mesh = model.getFEEngine().getMesh(); const auto & connectivity = mesh.getConnectivity(type); auto & displacement = model.getDisplacement(); Array update(displacement.size()); update.clear(); auto conn_it = connectivity.begin(connectivity.getNbComponent()); auto conn_end = connectivity.begin(connectivity.getNbComponent()); for (;conn_it != conn_end;++conn_it) { const auto & conn = *conn_it; for (UInt n = 0; n < conn.size(); ++n) { UInt node = conn(n); if (!update(node)) { Vector node_disp(displacement.storage() + node * spatial_dimension, spatial_dimension); node_disp += increment; update(node) = true; } } } } /* -------------------------------------------------------------------------- */ Checker::Checker(const SolidMechanicsModelCohesive & model, const Array & elements, ElementType type) - : model(model), type(std::move(type)), elements(elements), + : model(model), type(std::move(type)), //elements(elements), mat_cohesive(model.getMaterial(1)), sigma_c(mat_cohesive.get("sigma_c")), beta(mat_cohesive.get("beta")), G_c(mat_cohesive.get("G_c")), delta_0(mat_cohesive.get("delta_0")), kappa(mat_cohesive.get("kappa")), spatial_dimension(model.getSpatialDimension()), type_facet(Mesh::getFacetType(type)), type_cohesive(FEEngine::getCohesiveElementType(type_facet)), traction(mat_cohesive.getArray("tractions", type_cohesive)), damage(mat_cohesive.getArray("damage", type_cohesive)), nb_quad_per_el(model.getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(type_cohesive)), nb_element(model.getMesh().getNbElement(type_cohesive)), beta2_kappa2(beta * beta / kappa / kappa), beta2_kappa(beta * beta / kappa) { const Mesh & mesh = model.getMesh(); const auto & connectivity = mesh.getConnectivity(type); const auto & position = mesh.getNodes(); auto conn_it = connectivity.begin(connectivity.getNbComponent()); for (const auto & element : elements) { Vector conn_el(conn_it[element]); for (UInt n = 0; n < conn_el.size(); ++n) { UInt node = conn_el(n); if (Math::are_float_equal(position(node, _x), 0.)) nodes_to_check.insert(node); } } delta_c = 2 * G_c / sigma_c; sigma_c *= delta_c / (delta_c - delta_0); } /* -------------------------------------------------------------------------- */ void Checker::checkTractions(const Vector & opening, const Matrix & rotation) { auto normal_opening = opening * Vector{1., 0., 0.}; auto tangential_opening = opening - normal_opening; const Real normal_opening_norm = normal_opening.norm(); const Real tangential_opening_norm = tangential_opening.norm(); const Real delta = std::max(std::sqrt(tangential_opening_norm * tangential_opening_norm * beta2_kappa2 + normal_opening_norm * normal_opening_norm), 0.); Real theoretical_damage = std::min(delta / delta_c, 1.); theoretical_traction = (tangential_opening * beta2_kappa + normal_opening) * sigma_c / delta * (1. - theoretical_damage); // adjust damage theoretical_damage = std::max((delta - delta_0) / (delta_c - delta_0), 0.); theoretical_damage = std::min(theoretical_damage, 1.); Vector theoretical_traction_rotated = rotation * theoretical_traction; std::for_each( traction.begin(spatial_dimension), traction.end(spatial_dimension), [&theoretical_traction_rotated](auto && traction) { Real diff = Vector(theoretical_traction_rotated - traction).norm(); if (diff > 1e-14) throw std::domain_error("Tractions are incorrect"); }); std::for_each(damage.begin(), damage.end(), [&theoretical_damage](auto && damage) { if (not Math::are_float_equal(theoretical_damage, damage)) throw std::domain_error("Damage is incorrect"); }); } /* -------------------------------------------------------------------------- */ void Checker::computeEnergy(const Vector & opening) { /// compute energy auto Do = opening - opening_old; auto Dt = traction_old + theoretical_traction; Ed += .5 * Do.dot(Dt); opening_old = opening; traction_old = theoretical_traction; } /* -------------------------------------------------------------------------- */ void Checker::checkEquilibrium() { Vector residual_sum(spatial_dimension, 0.); const auto & residual = model.getInternalForce(); auto res_it = residual.begin(spatial_dimension); auto res_end = residual.end(spatial_dimension); for (; res_it != res_end; ++res_it) residual_sum += *res_it; if (!Math::are_float_equal(residual_sum.norm(), 0.)) throw std::domain_error("System is not in equilibrium!"); } /* -------------------------------------------------------------------------- */ void Checker::checkResidual(const Matrix & rotation) { Vector total_force(spatial_dimension, 0.); const auto & residual = model.getInternalForce(); for (auto node : nodes_to_check) { Vector res(residual.begin(spatial_dimension)[node]); total_force += res; } Vector theoretical_total_force(spatial_dimension); theoretical_total_force.mul(rotation, theoretical_traction); theoretical_total_force *= -1 * 2 * 2; for (UInt s = 0; s < spatial_dimension; ++s) { if (!Math::are_float_equal(total_force(s), theoretical_total_force(s))) { std::cout << "Total force isn't correct!" << std::endl; std::terminate(); } } }