diff --git a/src/common/aka_common.cc b/src/common/aka_common.cc index cc1df0c57..95b6c0598 100644 --- a/src/common/aka_common.cc +++ b/src/common/aka_common.cc @@ -1,151 +1,153 @@ /** * @file aka_common.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Wed Jan 13 2016 * * @brief Initialization of global variables * * @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 "aka_static_memory.hh" #include "static_communicator.hh" #include "static_solver.hh" #include "aka_random_generator.hh" #include "parser.hh" #include "cppargparse.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void initialize(int & argc, char **& argv) { AKANTU_DEBUG_IN(); initialize("", argc, argv); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void initialize(const std::string & input_file, int & argc, char **& argv) { AKANTU_DEBUG_IN(); StaticMemory::getStaticMemory(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(argc, argv); debug::debugger.setParallelContext(comm.whoAmI(), comm.getNbProc()); debug::initSignalHandler(); static_argparser.setParallelContext(comm.whoAmI(), comm.getNbProc()); static_argparser.setExternalExitFunction(debug::exit); static_argparser.addArgument("--aka_input_file", "Akantu's input file", 1, cppargparse::_string, std::string()); static_argparser.addArgument( "--aka_debug_level", std::string("Akantu's overall debug level") + std::string(" (0: error, 1: exceptions, 4: warnings, 5: info, ..., " "100: dump,") + std::string(" more info on levels can be foind in aka_error.hh)"), 1, cppargparse::_integer, int(dblWarning)); static_argparser.addArgument( "--aka_print_backtrace", "Should Akantu print a backtrace in case of error", 0, cppargparse::_boolean, false, true); static_argparser.parse(argc, argv, cppargparse::_remove_parsed); std::string infile = static_argparser["aka_input_file"]; if (infile == "") infile = input_file; debug::setDebugLevel(dblError); if ("" != infile) { readInputFile(infile); } long int seed; try { seed = static_parser.getParameter("seed", _ppsc_current_scope); } catch (debug::Exception &) { seed = time(NULL); } int dbl_level = static_argparser["aka_debug_level"]; debug::setDebugLevel(DebugLevel(dbl_level)); debug::debugger.printBacktrace(static_argparser["aka_print_backtrace"]); seed *= (comm.whoAmI() + 1); #if not defined(_WIN32) Rand48Generator::seed(seed); #endif RandGenerator::seed(seed); AKANTU_DEBUG_INFO("Random seed set to " << seed); /// initialize external solvers StaticSolver::getStaticSolver().initialize(argc, argv); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void finalize() { AKANTU_DEBUG_IN(); - if (StaticMemory::isInstantiated()) - delete &(StaticMemory::getStaticMemory()); if (StaticCommunicator::isInstantiated()) { StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); delete &comm; } + if (StaticMemory::isInstantiated()) { + delete &(StaticMemory::getStaticMemory()); + } + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void readInputFile(const std::string & input_file) { static_parser.parse(input_file); } /* -------------------------------------------------------------------------- */ cppargparse::ArgumentParser & getStaticArgumentParser() { return static_argparser; } /* -------------------------------------------------------------------------- */ Parser & getStaticParser() { return static_parser; } /* -------------------------------------------------------------------------- */ const ParserSection & getUserParser() { return *(static_parser.getSubSections(_st_user).first); } __END_AKANTU__ diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc index b35c54446..8bc2ed256 100644 --- a/src/common/aka_error.cc +++ b/src/common/aka_error.cc @@ -1,355 +1,355 @@ /** * @file aka_error.cc * * @author Nicolas Richart * * @date creation: Mon Sep 06 2010 * @date last modification: Mon Nov 16 2015 * * @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_config.hh" #include "aka_error.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #include #include -#if (defined(READLINK_COMMAND) || \ - defined(ADDR2LINE_COMMAND)) && \ +#if (defined(READLINK_COMMAND) || defined(ADDR2LINE_COMMAND)) && \ (not defined(_WIN32)) -# include -# include +#include +#include #endif #include #include #include #include #include #include #include #include #if defined(AKANTU_USE_OBSOLETE_GETTIMEOFDAY) -# include +#include #else -# include +#include #endif #ifdef AKANTU_USE_MPI -# include +#include #endif -#define __BEGIN_AKANTU_DEBUG__ namespace akantu { namespace debug { -#define __END_AKANTU_DEBUG__ } } +#define __BEGIN_AKANTU_DEBUG__ \ + namespace akantu { \ + namespace debug { +#define __END_AKANTU_DEBUG__ \ + } \ + } /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU_DEBUG__ static void printBacktraceAndExit(int sig) { printBacktrace(sig); debugger.exit(-50); } /* ------------------------------------------------------------------------ */ void initSignalHandler() { #if not defined(_WIN32) struct sigaction action; action.sa_handler = &printBacktraceAndExit; sigemptyset(&(action.sa_mask)); action.sa_flags = SA_RESETHAND; sigaction(SIGSEGV, &action, NULL); sigaction(SIGABRT, &action, NULL); #else std::signal(SIGSEGV, &printBacktraceAndExit); std::signal(SIGABRT, &printBacktraceAndExit); #endif } /* ------------------------------------------------------------------------ */ -std::string demangle(const char *symbol) { +std::string demangle(const char * symbol) { int status; std::string result; - char *demangled_name; + char * demangled_name; - if ((demangled_name = abi::__cxa_demangle(symbol, NULL, 0, &status)) != NULL) { + if ((demangled_name = abi::__cxa_demangle(symbol, NULL, 0, &status)) != + NULL) { result = demangled_name; free(demangled_name); } else { result = symbol; } return result; } /* ------------------------------------------------------------------------ */ -#if (defined(READLINK_COMMAND) || defined(ADDR2LINK_COMMAND)) && (not defined(_WIN32)) +#if (defined(READLINK_COMMAND) || defined(ADDR2LINK_COMMAND)) && \ + (not defined(_WIN32)) std::string exec(std::string cmd) { - FILE *pipe = popen(cmd.c_str(), "r"); - if (!pipe) return ""; + FILE * pipe = popen(cmd.c_str(), "r"); + if (!pipe) + return ""; char buffer[1024]; std::string result = ""; while (!feof(pipe)) { if (fgets(buffer, 128, pipe) != NULL) 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::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; + 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; + if (me != "") + addr_map[me] = 0; #endif - /// \todo for windows this part could be coded using CaptureStackBackTrace and SymFromAddr + /// \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; + 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; size_t 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) { + 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); + location = exec(std::string(BOOST_PP_STRINGIZE(READLINK_COMMAND)) + + std::string(" -f ") + location); #endif - std::string call = demangle(bt_line.substr(first + 1, second - first - 1).c_str()); + 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; + size_t addr; + sstra >> std::hex >> addr; std::cerr << location << " [" << call << "]"; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) - std::map ::iterator it = addr_map.find(location); + std::map::iterator 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; + 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; + 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 } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ Debugger::Debugger() { cout = &std::cerr; level = dblWarning; parallel_context = ""; file_open = false; print_backtrace = false; } /* ------------------------------------------------------------------------ */ Debugger::~Debugger() { if (file_open) { - dynamic_cast (cout)->close(); + dynamic_cast(cout)->close(); delete cout; } } - /* ------------------------------------------------------------------------ */ void Debugger::exit(int status) { if (status != EXIT_SUCCESS && status != -50) { #ifndef AKANTU_NDEBUG int * a = NULL; *a = 1; #endif - if(this->print_backtrace) + if (this->print_backtrace) akantu::debug::printBacktrace(15); } #ifdef AKANTU_USE_MPI if (status != EXIT_SUCCESS) MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN); #endif - std::exit(status); // not called when compiled with MPI due to MPI_Abort, but + std::exit( + status); // not called when compiled with MPI due to MPI_Abort, but // MPI_Abort does not have the noreturn attribute } /*------------------------------------------------------------------------- */ 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 - throw(akantu::debug::Exception) { + const std::string & file, unsigned int line, + __attribute__((unused)) bool silent, + __attribute__((unused)) + const std::string & location) const + throw(akantu::debug::Exception) { #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 { + const DebugLevel & level, + const std::string & info) const { if (this->level >= level) { #if 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 - << "{" << (unsigned int)timestamp << "} " - << prefix << " " << info - << std::endl; + *(cout) << parallel_context << "{" << (unsigned int)timestamp << "} " + << prefix << " " << info << std::endl; } } /* ------------------------------------------------------------------------ */ -void Debugger::setDebugLevel(const DebugLevel & level) { - this->level = level; -} +void Debugger::setDebugLevel(const DebugLevel & level) { this->level = level; } /* ------------------------------------------------------------------------ */ -const DebugLevel & Debugger::getDebugLevel() const { - return level; -} +const DebugLevel & Debugger::getDebugLevel() const { return level; } /* ------------------------------------------------------------------------ */ void Debugger::setLogFile(const std::string & filename) { if (file_open) { - dynamic_cast (cout)->close(); + dynamic_cast(cout)->close(); delete cout; } - std::ofstream *fileout = new std::ofstream(filename.c_str()); + std::ofstream * fileout = new std::ofstream(filename.c_str()); file_open = true; cout = fileout; } -std::ostream & Debugger::getOutputStream() { - return *cout; -} +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 << "] "; + 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); -} +void setDebugLevel(const DebugLevel & level) { debugger.setDebugLevel(level); } -const DebugLevel & getDebugLevel() { - return debugger.getDebugLevel(); -} +const DebugLevel & getDebugLevel() { return debugger.getDebugLevel(); } /* -------------------------------------------------------------------------- */ -void exit(int status) { - debugger.exit(status); -} +void exit(int status) { debugger.exit(status); } __END_AKANTU_DEBUG__ diff --git a/src/io/dumper/dumper_compute.hh b/src/io/dumper/dumper_compute.hh index aa5bd5595..8bba884d6 100644 --- a/src/io/dumper/dumper_compute.hh +++ b/src/io/dumper/dumper_compute.hh @@ -1,272 +1,273 @@ /** * @file dumper_compute.hh * * @author Guillaume Anciaux * * @date creation: Tue Sep 02 2014 * @date last modification: Sun Nov 15 2015 * * @brief Field that map a function to another 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 . * */ #ifndef __AKANTU_DUMPER_COMPUTE_HH__ #define __AKANTU_DUMPER_COMPUTE_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "dumper_iohelper.hh" #include "dumper_type_traits.hh" #include "dumper_field.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ __BEGIN_AKANTU_DUMPER__ class ComputeFunctorInterface { public: virtual ~ComputeFunctorInterface(){}; virtual UInt getDim() = 0; virtual UInt getNbComponent(UInt old_nb_comp) = 0; }; /* -------------------------------------------------------------------------- */ template class ComputeFunctorOutput : public ComputeFunctorInterface { public: ComputeFunctorOutput(){}; virtual ~ComputeFunctorOutput(){}; }; /* -------------------------------------------------------------------------- */ -template +template class ComputeFunctor : public ComputeFunctorOutput { public: ComputeFunctor(){}; virtual ~ComputeFunctor(){}; virtual return_type func(const input_type & d, Element global_index) = 0; }; /* -------------------------------------------------------------------------- */ template class FieldCompute : public Field { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ public: typedef typename SubFieldCompute::iterator sub_iterator; typedef typename SubFieldCompute::types sub_types; typedef typename sub_types::return_type sub_return_type; typedef _return_type return_type; typedef typename sub_types::data_type data_type; - typedef TypeTraits > types; - + typedef TypeTraits > + types; class iterator { public: - iterator(const sub_iterator & it, ComputeFunctor & func) - : it(it), func(func) {} - - bool operator!=(const iterator & it)const { return it.it != this->it; } - iterator operator++() { ++this->it; return *this; } - - - UInt currentGlobalIndex(){ - return this->it.currentGlobalIndex(); + iterator(const sub_iterator & it, + ComputeFunctor & func) + : it(it), func(func) {} + + bool operator!=(const iterator & it) const { return it.it != this->it; } + iterator operator++() { + ++this->it; + return *this; } - return_type operator*() { - return func.func(*it,it.getCurrentElement()); - } + UInt currentGlobalIndex() { return this->it.currentGlobalIndex(); } - Element getCurrentElement(){ - return this->it.getCurrentElement(); - } + return_type operator*() { return func.func(*it, it.getCurrentElement()); } + + Element getCurrentElement() { return this->it.getCurrentElement(); } UInt element_type() { return this->it.element_type(); } protected: sub_iterator it; - ComputeFunctor & func; + ComputeFunctor & func; }; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FieldCompute(SubFieldCompute & cont, ComputeFunctorInterface & func) - :sub_field(cont),func(dynamic_cast &>(func)){ + : sub_field(cont), + func(dynamic_cast &>( + func)) { this->checkHomogeneity(); }; + ~FieldCompute() { + delete &(this->sub_field); + delete &(this->func); + } + virtual void registerToDumper(const std::string & id, - iohelper::Dumper & dumper) { + iohelper::Dumper & dumper) { dumper.addElemDataField(id, *this); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: iterator begin() { return iterator(sub_field.begin(), func); } - iterator end () { return iterator(sub_field.end(), func); } + iterator end() { return iterator(sub_field.end(), func); } - UInt getDim() { - return func.getDim(); - } + UInt getDim() { return func.getDim(); } UInt size() { throw; // return Functor::size(); return 0; } - virtual void checkHomogeneity(){this->homogeneous = true;}; + virtual void checkHomogeneity() { this->homogeneous = true; }; - iohelper::DataType getDataType() { return iohelper::getDataType(); } + iohelper::DataType getDataType() { + return iohelper::getDataType(); + } /// get the number of components of the hosted field - virtual ElementTypeMap getNbComponents(UInt dim = _all_dimensions, - GhostType ghost_type = _not_ghost, - ElementKind kind = _ek_not_defined){ + virtual ElementTypeMap + getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, + ElementKind kind = _ek_not_defined) { ElementTypeMap nb_components; const ElementTypeMap & old_nb_components = - this->sub_field.getNbComponents(dim,ghost_type,kind); - + this->sub_field.getNbComponents(dim, ghost_type, kind); - ElementTypeMap::type_iterator tit = old_nb_components.firstType(dim,ghost_type,kind); - ElementTypeMap::type_iterator end = old_nb_components.lastType(dim,ghost_type,kind); + ElementTypeMap::type_iterator tit = + old_nb_components.firstType(dim, ghost_type, kind); + ElementTypeMap::type_iterator end = + old_nb_components.lastType(dim, ghost_type, kind); - while (tit != end){ - UInt nb_comp = old_nb_components(*tit,ghost_type); - nb_components(*tit,ghost_type) = func.getNbComponent(nb_comp); + while (tit != end) { + UInt nb_comp = old_nb_components(*tit, ghost_type); + nb_components(*tit, ghost_type) = func.getNbComponent(nb_comp); ++tit; } return nb_components; }; /// for connection to a FieldCompute inline virtual Field * connect(FieldComputeProxy & proxy); /// for connection to a FieldCompute virtual ComputeFunctorInterface * connect(HomogenizerProxy & proxy); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: SubFieldCompute & sub_field; - ComputeFunctor & func; + ComputeFunctor & func; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class FieldComputeProxy { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - - FieldComputeProxy(ComputeFunctorInterface & func):func(func){}; + FieldComputeProxy(ComputeFunctorInterface & func) : func(func){}; inline static Field * createFieldCompute(Field * field, - ComputeFunctorInterface & func){ + ComputeFunctorInterface & func) { + /// that looks fishy an object passed as a ref and destroyed at their end of + /// the function FieldComputeProxy compute_proxy(func); return field->connect(compute_proxy); } - template - Field * connectToField(T * ptr){ - if (dynamic_cast > *>(&func)){ + template Field * connectToField(T * ptr) { + if (dynamic_cast > *>(&func)) { return this->connectToFunctor >(ptr); - } - else if (dynamic_cast > *>(&func)){ + } else if (dynamic_cast > *>(&func)) { return this->connectToFunctor >(ptr); } - else if (dynamic_cast > *>(&func)){ + else if (dynamic_cast > *>(&func)) { return this->connectToFunctor >(ptr); } - else if (dynamic_cast > *>(&func)){ + else if (dynamic_cast > *>(&func)) { return this->connectToFunctor >(ptr); } - else throw; + else + throw; } - template - Field * connectToFunctor(T * ptr){ - FieldCompute * functor_ptr = new FieldCompute(*ptr,func); + template Field * connectToFunctor(T * ptr) { + FieldCompute * functor_ptr = + new FieldCompute(*ptr, func); return functor_ptr; } - template - Field * connectToFunctor(FieldCompute, - return_type2> * ptr){ - throw; // return new FieldCompute(*ptr,func); + template + Field * + connectToFunctor(FieldCompute, + return_type2> * ptr) { + throw; // return new FieldCompute(*ptr,func); return NULL; } - template + template Field * connectToFunctor(FieldCompute< - FieldCompute< - FieldCompute< - FieldCompute, - return_type2>, - return_type3>, - return_type4> * ptr){ - throw; // return new FieldCompute(*ptr,func); - return NULL; - } + FieldCompute, + return_type2>, + return_type3>, + return_type4> * ptr) { + throw; // return new FieldCompute(*ptr,func); + return NULL; + } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: ComputeFunctorInterface & func; - }; /* -------------------------------------------------------------------------- */ /// for connection to a FieldCompute template -inline Field * FieldCompute -::connect(FieldComputeProxy & proxy){ +inline Field * +FieldCompute::connect(FieldComputeProxy & proxy) { return proxy.connectToField(this); } /* -------------------------------------------------------------------------- */ __END_AKANTU_DUMPER__ __END_AKANTU__ - #endif /* __AKANTU_DUMPER_COMPUTE_HH__ */ diff --git a/src/io/dumper/dumper_field.hh b/src/io/dumper/dumper_field.hh index 12db7c388..68b40f12b 100644 --- a/src/io/dumper/dumper_field.hh +++ b/src/io/dumper/dumper_field.hh @@ -1,124 +1,127 @@ /** * @file dumper_field.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Sep 02 2014 * @date last modification: Tue Dec 08 2015 * * @brief Common interface for fields * * @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 . * */ #ifndef __AKANTU_DUMPER_FIELD_HH__ #define __AKANTU_DUMPER_FIELD_HH__ /* -------------------------------------------------------------------------- */ #include "dumper_iohelper.hh" /* -------------------------------------------------------------------------- */ - __BEGIN_AKANTU__ __BEGIN_AKANTU_DUMPER__ /* -------------------------------------------------------------------------- */ class FieldComputeProxy; class FieldComputeBaseInterface; class ComputeFunctorInterface; class HomogenizerProxy; /* -------------------------------------------------------------------------- */ - /// Field interface class Field { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - Field(): homogeneous(false) {} - virtual ~Field() {}; + Field() : homogeneous(false) {} + virtual ~Field(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: #ifdef AKANTU_USE_IOHELPER /// register this to the provided dumper virtual void registerToDumper(const std::string & id, iohelper::Dumper & dumper) = 0; #endif /// set the number of data per item (used for elements fields at the moment) - virtual void setNbData(UInt nb_data){AKANTU_DEBUG_TO_IMPLEMENT();}; + virtual void setNbData(UInt nb_data) { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// set the number of data per elem (used for elements fields at the moment) - virtual void setNbDataPerElem(const ElementTypeMap & nb_data){AKANTU_DEBUG_TO_IMPLEMENT();}; + virtual void setNbDataPerElem(const ElementTypeMap & nb_data) { + AKANTU_DEBUG_TO_IMPLEMENT(); + }; /// set the number of data per elem (used for elements fields at the moment) - virtual void setNbDataPerElem(UInt nb_data){AKANTU_DEBUG_TO_IMPLEMENT();}; + virtual void setNbDataPerElem(UInt nb_data) { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// get the number of components of the hosted field - virtual ElementTypeMap getNbComponents(UInt dim = _all_dimensions, - GhostType ghost_type = _not_ghost, - ElementKind kind = _ek_not_defined){throw;}; + virtual ElementTypeMap + getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, + ElementKind kind = _ek_not_defined) { + throw; + }; /// for connection to a FieldCompute - inline virtual Field * connect(FieldComputeProxy & proxy){throw;}; + inline virtual Field * connect(FieldComputeProxy & proxy) { throw; }; /// for connection to a FieldCompute - inline virtual ComputeFunctorInterface * connect(HomogenizerProxy & proxy){throw;}; + inline virtual ComputeFunctorInterface * connect(HomogenizerProxy & proxy) { + throw; + }; /// check if the same quantity of data for all element types virtual void checkHomogeneity() = 0; /// return the dumper name - std::string getGroupName(){return group_name;}; + std::string getGroupName() { return group_name; }; /// return the id of the field - std::string getID(){return field_id;}; + std::string getID() { return field_id; }; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the flag to know if the field is homogeneous/contiguous virtual bool isHomogeneous() { return homogeneous; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the flag to know if it is homogeneous bool homogeneous; /// the name of the group it was associated to std::string group_name; /// the name of the dumper it was associated to std::string field_id; }; /* -------------------------------------------------------------------------- */ __END_AKANTU_DUMPER__ __END_AKANTU__ - #endif /* __AKANTU_DUMPER_FIELD_HH__ */ diff --git a/src/io/dumper/dumper_generic_elemental_field.hh b/src/io/dumper/dumper_generic_elemental_field.hh index 52074e644..0302a9022 100644 --- a/src/io/dumper/dumper_generic_elemental_field.hh +++ b/src/io/dumper/dumper_generic_elemental_field.hh @@ -1,224 +1,224 @@ /** * @file dumper_generic_elemental_field.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Sep 02 2014 * @date last modification: Mon Sep 21 2015 * * @brief Generic interface for elemental fields * * @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 . * */ #ifndef __AKANTU_DUMPER_GENERIC_ELEMENTAL_FIELD_HH__ #define __AKANTU_DUMPER_GENERIC_ELEMENTAL_FIELD_HH__ /* -------------------------------------------------------------------------- */ #include "dumper_field.hh" #include "element_type_map_filter.hh" #include "dumper_element_iterator.hh" #include "dumper_homogenizing_field.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ __BEGIN_AKANTU_DUMPER__ /* -------------------------------------------------------------------------- */ template class iterator_type> class GenericElementalField : public Field { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ public: // check dumper_type_traits.hh for additional information over these types typedef _types types; typedef typename types::data_type data_type; - typedef typename types::it_type it_type; + typedef typename types::it_type it_type; typedef typename types::field_type field_type; typedef typename types::array_type array_type; typedef typename types::array_iterator array_iterator; typedef typename field_type::type_iterator field_type_iterator; typedef iterator_type iterator; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GenericElementalField(const field_type & field, UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) : field(field), spatial_dimension(spatial_dimension), ghost_type(ghost_type), element_kind(element_kind) { this->checkHomogeneity(); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get the number of components of the hosted field virtual ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) { return this->field.getNbComponents(dim, ghost_type, kind); }; /// return the size of the contained data: i.e. the number of elements ? virtual UInt size() { checkHomogeneity(); return this->nb_total_element; } /// return the iohelper datatype to be dumped iohelper::DataType getDataType() { return iohelper::getDataType(); } protected: /// return the number of entries per element UInt getNbDataPerElem(const ElementType & type, const GhostType & ghost_type = _not_ghost) const { if (!nb_data_per_elem.exists(type, ghost_type)) return field(type, ghost_type).getNbComponent(); return nb_data_per_elem(type, this->ghost_type); } /// check if the same quantity of data for all element types virtual void checkHomogeneity(); public: virtual void registerToDumper(const std::string & id, iohelper::Dumper & dumper) { dumper.addElemDataField(id, *this); }; /// for connection to a FieldCompute inline virtual Field * connect(FieldComputeProxy & proxy) { return proxy.connectToField(this); } /// for connection to a Homogenizer inline virtual ComputeFunctorInterface * connect(HomogenizerProxy & proxy) { return proxy.connectToField(this); }; virtual iterator begin() { field_type_iterator tit; field_type_iterator end; /// type iterators on the elemental field tit = this->field.firstType(this->spatial_dimension, this->ghost_type, this->element_kind); end = this->field.lastType(this->spatial_dimension, this->ghost_type, this->element_kind); /// skip all types without data ElementType type = *tit; for (; tit != end && this->field(*tit, this->ghost_type).getSize() == 0; ++tit) { } type = *tit; if (tit == end) return this->end(); /// getting information for the field of the given type const array_type & vect = this->field(type, this->ghost_type); UInt nb_data_per_elem = this->getNbDataPerElem(type); UInt nb_component = vect.getNbComponent(); UInt size = (vect.getSize() * nb_component) / nb_data_per_elem; /// define element-wise iterator array_iterator it = vect.begin_reinterpret(nb_data_per_elem, size); array_iterator it_end = vect.end_reinterpret(nb_data_per_elem, size); /// define data iterator iterator rit = iterator(this->field, tit, end, it, it_end, this->ghost_type); rit.setNbDataPerElem(this->nb_data_per_elem); return rit; } virtual iterator end() { field_type_iterator tit; field_type_iterator end; tit = this->field.firstType(this->spatial_dimension, this->ghost_type, this->element_kind); end = this->field.lastType(this->spatial_dimension, this->ghost_type, this->element_kind); ElementType type = *tit; for (; tit != end; ++tit) type = *tit; const array_type & vect = this->field(type, this->ghost_type); UInt nb_data = this->getNbDataPerElem(type); UInt nb_component = vect.getNbComponent(); UInt size = (vect.getSize() * nb_component) / nb_data; array_iterator it = vect.end_reinterpret(nb_data, size); iterator rit = iterator(this->field, end, end, it, it, this->ghost_type); rit.setNbDataPerElem(this->nb_data_per_elem); return rit; } virtual UInt getDim() { if (this->homogeneous) { field_type_iterator tit = this->field.firstType( this->spatial_dimension, this->ghost_type, this->element_kind); return this->getNbDataPerElem(*tit); } throw; return 0; } void setNbDataPerElem(const ElementTypeMap & nb_data) { nb_data_per_elem = nb_data; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the ElementTypeMapArray embedded in the field const field_type & field; /// total number of elements UInt nb_total_element; /// the spatial dimension of the problem UInt spatial_dimension; /// whether this is a ghost field or not (for type selection) GhostType ghost_type; /// The element kind to operate on ElementKind element_kind; /// The number of data per element type ElementTypeMap nb_data_per_elem; }; /* -------------------------------------------------------------------------- */ #include "dumper_generic_elemental_field_tmpl.hh" /* -------------------------------------------------------------------------- */ __END_AKANTU_DUMPER__ __END_AKANTU__ #endif /* __AKANTU_DUMPER_GENERIC_ELEMENTAL_FIELD_HH__ */ diff --git a/src/io/dumper/dumper_iohelper.hh b/src/io/dumper/dumper_iohelper.hh index c905a9925..72668420b 100644 --- a/src/io/dumper/dumper_iohelper.hh +++ b/src/io/dumper/dumper_iohelper.hh @@ -1,157 +1,158 @@ /** * @file dumper_iohelper.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Oct 26 2012 * @date last modification: Fri Apr 17 2015 * * @brief Define the akantu dumper interface for IOhelper dumpers * * @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 "aka_types.hh" #include "aka_array.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DUMPER_IOHELPER_HH__ #define __AKANTU_DUMPER_IOHELPER_HH__ /* -------------------------------------------------------------------------- */ - namespace iohelper { - class Dumper; +class Dumper; } __BEGIN_AKANTU__ UInt getIOHelperType(ElementType type); namespace dumper { - class Field; - class VariableBase; +class Field; +class VariableBase; } class Mesh; class DumperIOHelper { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DumperIOHelper(); virtual ~DumperIOHelper(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register a given Mesh for the current dumper - virtual void registerMesh(const Mesh & mesh, UInt spatial_dimension = _all_dimensions, - const GhostType & ghost_type = _not_ghost, - const ElementKind & element_kind = _ek_not_defined); + virtual void registerMesh(const Mesh & mesh, + UInt spatial_dimension = _all_dimensions, + const GhostType & ghost_type = _not_ghost, + const ElementKind & element_kind = _ek_not_defined); /// register a filtered Mesh (provided filter lists) for the current dumper - virtual void registerFilteredMesh(const Mesh & mesh, - const ElementTypeMapArray & elements_filter, - const Array & nodes_filter, - UInt spatial_dimension = _all_dimensions, - const GhostType & ghost_type = _not_ghost, - const ElementKind & element_kind = _ek_not_defined); + virtual void + registerFilteredMesh(const Mesh & mesh, + const ElementTypeMapArray & elements_filter, + const Array & nodes_filter, + UInt spatial_dimension = _all_dimensions, + const GhostType & ghost_type = _not_ghost, + const ElementKind & element_kind = _ek_not_defined); /// register a Field object identified by name and provided by pointer void registerField(const std::string & field_id, dumper::Field * field); /// remove the Field identified by name from managed fields void unRegisterField(const std::string & field_id); /// register a VariableBase object identified by name and provided by pointer - void registerVariable(const std::string & variable_id, dumper::VariableBase * variable); + void registerVariable(const std::string & variable_id, + dumper::VariableBase * variable); /// remove a VariableBase identified by name from managed fields void unRegisterVariable(const std::string & variable_id); - /// request dump: this calls IOHelper dump routine + /// request dump: this calls IOHelper dump routine virtual void dump(); - /// request dump: this first set the current step and then calls IOHelper dump routine + /// request dump: this first set the current step and then calls IOHelper dump + /// routine virtual void dump(UInt step); - /// request dump: this first set the current step and current time and then calls IOHelper dump routine + /// request dump: this first set the current step and current time and then + /// calls IOHelper dump routine virtual void dump(Real current_time, UInt step); /// set the parallel context for IOHeper virtual void setParallelContext(bool is_parallel); /// set the directory where to generate the dumped files virtual void setDirectory(const std::string & directory); /// set the base name (needed by most IOHelper dumpers) virtual void setBaseName(const std::string & basename); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - /// direct access to the iohelper::Dumper object AKANTU_GET_MACRO(Dumper, *dumper, iohelper::Dumper &) /// set the timestep of the iohelper::Dumper void setTimeStep(Real time_step); public: - /* ------------------------------------------------------------------------ */ /* Variable wrapper */ - template::value> - class Variable; + template ::value> class Variable; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// internal iohelper::Dumper iohelper::Dumper * dumper; typedef std::map Fields; typedef std::map Variables; /// list of registered fields to dump Fields fields; Variables variables; /// dump counter UInt count; /// directory name std::string directory; /// filename prefix std::string filename; /// is time tracking activated in the dumper bool time_activated; }; __END_AKANTU__ #endif /* __AKANTU_DUMPER_IOHELPER_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index f27300a57..8ed212402 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1836 +1,1909 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Clement Roux * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Mon Nov 30 2015 * * @brief Implementation of the SolidMechanicsModel 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 "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "group_manager_inline_impl.cc" #include "dumpable_inline_impl.hh" #include "integration_scheme_2nd_order.hh" #include "element_group.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_PETSC #include "solver_petsc.hh" #include "petsc_matrix.hh" #endif #ifdef AKANTU_USE_IOHELPER -# include "dumper_field.hh" -# include "dumper_paraview.hh" -# include "dumper_homogenizing_field.hh" -# include "dumper_internal_material_field.hh" -# include "dumper_elemental_field.hh" -# include "dumper_material_padders.hh" -# include "dumper_element_partition.hh" -# include "dumper_iohelper.hh" +#include "dumper_field.hh" +#include "dumper_paraview.hh" +#include "dumper_homogenizing_field.hh" +#include "dumper_internal_material_field.hh" +#include "dumper_elemental_field.hh" +#include "dumper_material_padders.hh" +#include "dumper_element_partition.hh" +#include "dumper_iohelper.hh" #endif #ifdef AKANTU_DAMAGE_NON_LOCAL -# include "non_local_manager.hh" +#include "non_local_manager.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ -const SolidMechanicsModelOptions default_solid_mechanics_model_options(_explicit_lumped_mass, false); +const SolidMechanicsModelOptions + default_solid_mechanics_model_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ -SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, - UInt dim, - const ID & id, - const MemoryID & memory_id) : - Model(mesh, dim, id, memory_id), - BoundaryCondition(), - time_step(NAN), f_m2a(1.0), - mass_matrix(NULL), - velocity_damping_matrix(NULL), - stiffness_matrix(NULL), - jacobian_matrix(NULL), - material_index("material index", id), - material_local_numbering("material local numbering", id), - material_selector(new DefaultMaterialSelector(material_index)), - is_default_material_selector(true), - integrator(NULL), - increment_flag(false), solver(NULL), - synch_parallel(NULL), - are_materials_instantiated(false), - non_local_manager(NULL), - pbc_synch(NULL) { +SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, + const MemoryID & memory_id) + : Model(mesh, dim, id, memory_id), BoundaryCondition(), + time_step(NAN), f_m2a(1.0), mass_matrix(NULL), + velocity_damping_matrix(NULL), stiffness_matrix(NULL), + jacobian_matrix(NULL), material_index("material index", id), + material_local_numbering("material local numbering", id), + material_selector(new DefaultMaterialSelector(material_index)), + is_default_material_selector(true), integrator(NULL), + increment_flag(false), solver(NULL), synch_parallel(NULL), + are_materials_instantiated(false), non_local_manager(NULL), + pbc_synch(NULL) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); - registerFEEngineObject("SolidMechanicsFEEngine", mesh, spatial_dimension); + registerFEEngineObject("SolidMechanicsFEEngine", mesh, + spatial_dimension); this->displacement = NULL; - this->mass = NULL; - this->velocity = NULL; + this->mass = NULL; + this->velocity = NULL; this->acceleration = NULL; - this->force = NULL; - this->residual = NULL; + this->force = NULL; + this->residual = NULL; this->blocked_dofs = NULL; - this->increment = NULL; + this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; this->previous_displacement = NULL; materials.clear(); mesh.registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("paraview_all", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif AKANTU_DEBUG_OUT(); } - - /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; delete solver; delete mass_matrix; delete velocity_damping_matrix; - if(stiffness_matrix && stiffness_matrix != jacobian_matrix) + if (stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; delete jacobian_matrix; delete synch_parallel; - if(is_default_material_selector) { + if (is_default_material_selector) { delete material_selector; material_selector = NULL; } + flatten_internal_map::iterator fl_it = this->registered_internals.begin(); + flatten_internal_map::iterator fl_end = this->registered_internals.end(); + for (; fl_it != fl_end; ++fl_it) { + delete fl_it->second; + } + #ifdef AKANTU_DAMAGE_NON_LOCAL delete non_local_manager; non_local_manager = NULL; #endif delete pbc_synch; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(const ModelOptions & options) { Model::initFull(options); const SolidMechanicsModelOptions & smm_options = - dynamic_cast(options); + dynamic_cast(options); this->method = smm_options.analysis_method; // initialize the vectors this->initArrays(); // set the initial condition to 0 this->force->clear(); this->velocity->clear(); this->acceleration->clear(); this->displacement->clear(); // initialize pbc - if(this->pbc_pair.size()!=0) + if (this->pbc_pair.size() != 0) this->initPBC(); // initialize the time integration schemes - switch(this->method) { + switch (this->method) { case _explicit_lumped_mass: this->initExplicit(); break; case _explicit_consistent_mass: this->initSolver(); this->initExplicit(); break; case _implicit_dynamic: this->initImplicit(true); break; case _static: this->initImplicit(false); this->initArraysPreviousDisplacment(); break; default: AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel"); break; } #ifdef AKANTU_DAMAGE_NON_LOCAL /// create the non-local manager object for non-local damage computations this->non_local_manager = new NonLocalManager(*this); #endif // initialize the materials - if(this->parser->getLastParsedFile() != "") { + if (this->parser->getLastParsedFile() != "") { this->instantiateMaterials(); } - if(!smm_options.no_init_materials) { + if (!smm_options.no_init_materials) { this->initMaterials(); } - if(this->increment_flag) + if (this->increment_flag) this->initBC(*this, *this->displacement, *this->increment, *this->force); else this->initBC(*this, *this->displacement, *this->force); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, - DataAccessor * data_accessor) { + DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); - if (data_accessor == NULL) data_accessor = this; - synch_parallel = &createParallelSynch(partition,data_accessor); + if (data_accessor == NULL) + data_accessor = this; + synch_parallel = &createParallelSynch(partition, data_accessor); synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); synch_registry->registerSynchronizer(*synch_parallel, _gst_for_dump); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } - /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) { AKANTU_DEBUG_IN(); - //in case of switch from implicit to explicit - if(!this->isExplicit()) + // in case of switch from implicit to explicit + if (!this->isExplicit()) method = analysis_method; - if (integrator) delete integrator; + if (integrator) + delete integrator; integrator = new CentralDifference(); UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); - std::stringstream sstr; sstr << id << ":increment_acceleration"; - increment_acceleration = &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); + std::stringstream sstr; + sstr << id << ":increment_acceleration"; + increment_acceleration = + &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initArraysPreviousDisplacment() { AKANTU_DEBUG_IN(); this->setIncrementFlagOn(); - if(not this->previous_displacement) { + if (not this->previous_displacement) { UInt nb_nodes = this->mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << this->id << ":previous_displacement"; - this->previous_displacement = &(this->alloc (sstr_disp_t.str(), nb_nodes, this->spatial_dimension, 0.)); + this->previous_displacement = &(this->alloc( + sstr_disp_t.str(), nb_nodes, this->spatial_dimension, 0.)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); - std::stringstream sstr_disp; sstr_disp << id << ":displacement"; + std::stringstream sstr_disp; + sstr_disp << id << ":displacement"; // std::stringstream sstr_mass; sstr_mass << id << ":mass"; - std::stringstream sstr_velo; sstr_velo << id << ":velocity"; - std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; - std::stringstream sstr_forc; sstr_forc << id << ":force"; - std::stringstream sstr_resi; sstr_resi << id << ":residual"; - std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs"; - - displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); - // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); - velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); - acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); - force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); - residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); - blocked_dofs = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); - - std::stringstream sstr_curp; sstr_curp << id << ":current_position"; - current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); - - for(UInt g = _not_ghost; g <= _ghost; ++g) { - GhostType gt = (GhostType) g; - Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); - Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); - for(; it != end; ++it) { + std::stringstream sstr_velo; + sstr_velo << id << ":velocity"; + std::stringstream sstr_acce; + sstr_acce << id << ":acceleration"; + std::stringstream sstr_forc; + sstr_forc << id << ":force"; + std::stringstream sstr_resi; + sstr_resi << id << ":residual"; + std::stringstream sstr_boun; + sstr_boun << id << ":blocked_dofs"; + + displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, + REAL_INIT_VALUE)); + // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, + // 0)); + velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, + REAL_INIT_VALUE)); + acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, + REAL_INIT_VALUE)); + force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, + REAL_INIT_VALUE)); + residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, + REAL_INIT_VALUE)); + blocked_dofs = + &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); + + std::stringstream sstr_curp; + sstr_curp << id << ":current_position"; + current_position = + &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); + + for (UInt g = _not_ghost; g <= _ghost; ++g) { + GhostType gt = (GhostType)g; + Mesh::type_iterator it = + mesh.firstType(spatial_dimension, gt, _ek_not_defined); + Mesh::type_iterator end = + mesh.lastType(spatial_dimension, gt, _ek_not_defined); + for (; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); material_index.alloc(nb_element, 1, *it, gt); material_local_numbering.alloc(nb_element, 1, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initPBC() { Model::initPBC(); registerPBCSynchronizer(); - // as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves + // as long as there are ones on the diagonal of the matrix, we can put + // boudandary true for slaves std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); - while(it != end) { - for (UInt i=0; iregisterSynchronizer(*pbc_synch, _gst_smm_uv); synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_mass); synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_res); synch_registry->registerSynchronizer(*pbc_synch, _gst_for_dump); changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); Real * current_position_val = current_position->storage(); - Real * position_val = mesh.getNodes().storage(); - Real * displacement_val = displacement->storage(); + Real * position_val = mesh.getNodes().storage(); + Real * displacement_val = displacement->storage(); /// compute current_position = initial_position + displacement - memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); - for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { + memcpy(current_position_val, position_val, + nb_nodes * spatial_dimension * sizeof(Real)); + for (UInt n = 0; n < nb_nodes * spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions - memcpy(residual->storage(), force->storage(), nb_nodes*spatial_dimension*sizeof(Real)); + memcpy(residual->storage(), force->storage(), + nb_nodes * spatial_dimension * sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /*----------------------------------------------------------------------------*/ -void SolidMechanicsModel::reInitialize() -{ - -} - +void SolidMechanicsModel::reInitialize() {} /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * This function compute the second member of the motion equation. That is to * say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given * by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$ * F_{int} = \int_{\Omega} N \sigma d\Omega@f$ * */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); // f = f_ext - f_int // f = f_ext - if(need_initialize) initializeUpdateResidualData(); + if (need_initialize) + initializeUpdateResidualData(); AKANTU_DEBUG_INFO("Compute local stresses"); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ this->non_local_manager->computeAllNonLocalStresses(); #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the stress AKANTU_DEBUG_INFO("Send data for residual assembly"); synch_registry->asynchronousSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for local elements"); - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_not_ghost); } - AKANTU_DEBUG_INFO("Wait distant stresses"); // finalize communications synch_registry->waitEndSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeStresses() { if (isExplicit()) { // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); // compute stresses on all local elements for each materials std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } - /* ------------------------------------------------------------------------ */ +/* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ this->non_local_manager->computeAllNonLocalStresses(); #endif } else { std::vector::iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStressesFromTangentModuli(_not_ghost); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidualInternal() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Update the residual"); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; - if(method != _static) { + if (method != _static) { // f -= Ma - if(mass_matrix) { + if (mass_matrix) { // if full mass_matrix Array * Ma = new Array(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else if (mass) { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); - Real * mass_val = mass->storage(); - Real * accel_val = acceleration->storage(); - Real * res_val = residual->storage(); + Real * mass_val = mass->storage(); + Real * accel_val = acceleration->storage(); + Real * res_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { - if(!(*blocked_dofs_val)) { - *res_val -= *accel_val * *mass_val /f_m2a; - } - blocked_dofs_val++; - res_val++; - mass_val++; - accel_val++; + if (!(*blocked_dofs_val)) { + *res_val -= *accel_val ** mass_val / f_m2a; + } + blocked_dofs_val++; + res_val++; + mass_val++; + accel_val++; } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } // f -= Cv - if(velocity_damping_matrix) { + if (velocity_damping_matrix) { Array * Cv = new Array(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); updateResidualInternal(); - if(method == _explicit_lumped_mass) { + if (method == _explicit_lumped_mass) { /* residual = residual_{n+1} - M * acceleration_n therefore solution = increment acceleration not acceleration */ - solveLumped(*increment_acceleration, - *mass, - *residual, - *blocked_dofs, - f_m2a); + solveLumped(*increment_acceleration, *mass, *residual, *blocked_dofs, + f_m2a); } else if (method == _explicit_consistent_mass) { solve(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::solveLumped(Array & x, - const Array & A, - const Array & b, - const Array & blocked_dofs, - Real alpha) { +void SolidMechanicsModel::solveLumped(Array & x, const Array & A, + const Array & b, + const Array & blocked_dofs, + Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { - if(!(*blocked_dofs_val)) { - *x_val = alpha * (*b_val / *A_val); + if (!(*blocked_dofs_val)) { + *x_val = alpha *(*b_val / *A_val); } x_val++; A_val++; b_val++; blocked_dofs_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); - if(increment_flag) { - if(previous_displacement) increment->copy(*previous_displacement); - else increment->copy(*displacement); + if (increment_flag) { + if (previous_displacement) + increment->copy(*previous_displacement); + else + increment->copy(*displacement); } - AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " - << "have called initExplicit ? " - << "or initImplicit ?"); + AKANTU_DEBUG_ASSERT(integrator, "itegrator should have been allocated: " + << "have called initExplicit ? " + << "or initImplicit ?"); - integrator->integrationSchemePred(time_step, - *displacement, - *velocity, - *acceleration, - *blocked_dofs); + integrator->integrationSchemePred(time_step, *displacement, *velocity, + *acceleration, *blocked_dofs); - if(increment_flag) { + if (increment_flag) { Real * inc_val = increment->storage(); Real * dis_val = displacement->storage(); - UInt nb_degree_of_freedom = displacement->getSize() * displacement->getNbComponent(); + UInt nb_degree_of_freedom = + displacement->getSize() * displacement->getNbComponent(); for (UInt n = 0; n < nb_degree_of_freedom; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); - integrator->integrationSchemeCorrAccel(time_step, - *displacement, - *velocity, - *acceleration, - *blocked_dofs, - *increment_acceleration); + integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, + *acceleration, *blocked_dofs, + *increment_acceleration); - if(previous_displacement) previous_displacement->copy(*displacement); + if (previous_displacement) + previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStep() { AKANTU_DEBUG_IN(); - EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); + EventManager::sendEvent( + SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); this->explicitPred(); this->updateResidual(); this->updateAcceleration(); this->explicitCorr(); - EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method)); + EventManager::sendEvent( + SolidMechanicsModelEvent::AfterSolveStepEvent(method)); AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Initialize the solver and create the sparse matrices needed. * */ -void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { -#if !defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake +void SolidMechanicsModel::initSolver(__attribute__((unused)) + SolverOptions & options) { +#if !defined(AKANTU_USE_MUMPS) && \ + !defined(AKANTU_USE_PETSC) // or other solver in the future \todo add + // AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_nodes = mesh.getNbGlobalNodes(); delete jacobian_matrix; - std::stringstream sstr; sstr << id << ":jacobian_matrix"; + std::stringstream sstr; + sstr << id << ":jacobian_matrix"; #ifdef AKANTU_USE_PETSC - jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); + jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, + _symmetric, sstr.str(), memory_id); #else - jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _unsymmetric, sstr.str(), memory_id); -#endif //AKANTU_USE PETSC + jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, + _unsymmetric, sstr.str(), memory_id); +#endif // AKANTU_USE PETSC jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); if (!isExplicit()) { delete stiffness_matrix; - std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; + std::stringstream sstr_sti; + sstr_sti << id << ":stiffness_matrix"; #ifdef AKANTU_USE_PETSC - stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); + stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, + _symmetric, sstr.str(), memory_id); stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); #else - stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); -#endif //AKANTU_USE_PETSC + stiffness_matrix = + new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); +#endif // AKANTU_USE_PETSC } delete solver; - std::stringstream sstr_solv; sstr_solv << id << ":solver"; + std::stringstream sstr_solv; + sstr_solv << id << ":solver"; #ifdef AKANTU_USE_PETSC solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str()); #elif defined(AKANTU_USE_MUMPS) solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); -#endif //AKANTU_USE_MUMPS +#endif // AKANTU_USE_MUMPS - if(solver) + if (solver) solver->initialize(options); -#endif //AKANTU_HAS_SOLVER +#endif // AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initJacobianMatrix() { #if defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC) // @todo make it more flexible: this is an ugly patch to treat the case of non // fix profile of the K matrix delete jacobian_matrix; - std::stringstream sstr_sti; sstr_sti << id << ":jacobian_matrix"; - jacobian_matrix = new SparseMatrix(*stiffness_matrix, sstr_sti.str(), memory_id); + std::stringstream sstr_sti; + sstr_sti << id << ":jacobian_matrix"; + jacobian_matrix = + new SparseMatrix(*stiffness_matrix, sstr_sti.str(), memory_id); - std::stringstream sstr_solv; sstr_solv << id << ":solver"; + std::stringstream sstr_solv; + sstr_solv << id << ":solver"; delete solver; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); - if(solver) + if (solver) solver->initialize(_solver_no_options); #else AKANTU_DEBUG_ERROR("You need to activate the solver mumps."); #endif } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ -void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) { +void SolidMechanicsModel::initImplicit(bool dynamic, + SolverOptions & solver_options) { AKANTU_DEBUG_IN(); method = dynamic ? _implicit_dynamic : _static; - if (!increment) setIncrementFlagOn(); + if (!increment) + setIncrementFlagOn(); initSolver(solver_options); - if(method == _implicit_dynamic) { - if(integrator) delete integrator; + if (method == _implicit_dynamic) { + if (integrator) + delete integrator; integrator = new TrapezoidalRule2(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; - std::stringstream sstr; sstr << id << ":tmp_mass_matrix"; - SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id); + std::stringstream sstr; + sstr << id << ":tmp_mass_matrix"; + SparseMatrix * tmp_mass = + new SparseMatrix(*mass_matrix, sstr.str(), memory_id); #ifdef AKANTU_USE_MUMPS - std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix"; + std::stringstream sstr_solver; + sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); -#endif //AKANTU_USE_MUMPS +#endif // AKANTU_USE_MUMPS acc_solver->initialize(); tmp_mass->applyBoundary(*blocked_dofs); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); stiffness_matrix->clear(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { - if(!velocity_damping_matrix) - velocity_damping_matrix = - new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id); + if (!velocity_damping_matrix) + velocity_damping_matrix = new SparseMatrix( + *jacobian_matrix, id + ":velocity_damping_matrix", memory_id); return *velocity_damping_matrix; } /* -------------------------------------------------------------------------- */ -template<> -bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){ +template <> +bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, + Real & error) { AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); error = 0; Real norm[2] = {0., 0.}; - Real * increment_val = increment->storage(); - bool * blocked_dofs_val = blocked_dofs->storage(); + Real * increment_val = increment->storage(); + bool * blocked_dofs_val = blocked_dofs->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { - if(!(*blocked_dofs_val) && is_local_node) { - norm[0] += *increment_val * *increment_val; - norm[1] += *displacement_val * *displacement_val; + if (!(*blocked_dofs_val) && is_local_node) { + norm[0] += *increment_val * *increment_val; + norm[1] += *displacement_val * *displacement_val; } blocked_dofs_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); - AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); + AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), + "Something goes wrong in the solve phase"); if (norm[1] < Math::getTolerance()) { error = norm[0]; AKANTU_DEBUG_OUT(); // cout<<"Error 1: "< Math::getTolerance()) + if (norm[1] > Math::getTolerance()) error = norm[0] / norm[1]; else - error = norm[0]; //In case the total displacement is zero! + error = norm[0]; // In case the total displacement is zero! // cout<<"Error 2: "< -bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) { +template <> +bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, + Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); norm = 0; Real * residual_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); - if(is_local_node) { + if (is_local_node) { for (UInt d = 0; d < nb_degree_of_freedom; ++d) { - if(!(*blocked_dofs_val)) { - norm += *residual_val * *residual_val; - } - blocked_dofs_val++; - residual_val++; + if (!(*blocked_dofs_val)) { + norm += *residual_val * *residual_val; + } + blocked_dofs_val++; + residual_val++; } } else { blocked_dofs_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); - AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); + AKANTU_DEBUG_ASSERT(!Math::isnan(norm), + "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ -template<> -bool SolidMechanicsModel::testConvergence<_scc_residual_mass_wgh>(Real tolerance, - Real & norm) { +template <> +bool SolidMechanicsModel::testConvergence<_scc_residual_mass_wgh>( + Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); - - UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->storage(); Real * mass_val = this->mass->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); - if(is_local_node) { + if (is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { - if(!(*blocked_dofs_val)) { - norm += *residual_val * *residual_val/(*mass_val * *mass_val); - } - blocked_dofs_val++; - residual_val++; - mass_val++; + if (!(*blocked_dofs_val)) { + norm += *residual_val * *residual_val / (*mass_val * *mass_val); + } + blocked_dofs_val++; + residual_val++; + mass_val++; } } else { blocked_dofs_val += spatial_dimension; residual_val += spatial_dimension; mass_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); - AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); + AKANTU_DEBUG_ASSERT(!Math::isnan(norm), + "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ -bool SolidMechanicsModel::testConvergenceResidual(Real tolerance){ +bool SolidMechanicsModel::testConvergenceResidual(Real tolerance) { AKANTU_DEBUG_IN(); - Real error=0; + Real error = 0; bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ -bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error){ +bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, + Real & error) { AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ -bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance){ +bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance) { AKANTU_DEBUG_IN(); - Real error=0; + Real error = 0; bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ -bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error){ +bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, + Real & error) { AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); - if(method == _implicit_dynamic) - integrator->integrationSchemePred(time_step, - *displacement, - *velocity, - *acceleration, - *blocked_dofs); + if (method == _implicit_dynamic) + integrator->integrationSchemePred(time_step, *displacement, *velocity, + *acceleration, *blocked_dofs); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); - if(method == _implicit_dynamic) { - integrator->integrationSchemeCorrDispl(time_step, - *displacement, - *velocity, - *acceleration, - *blocked_dofs, - *increment); + if (method == _implicit_dynamic) { + integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, + *acceleration, *blocked_dofs, + *increment); } else { UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); bool * boun_val = blocked_dofs->storage(); - for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){ + for (UInt j = 0; j < nb_degree_of_freedom; + ++j, ++disp_val, ++incr_val, ++boun_val) { *incr_val *= (1. - *boun_val); *disp_val += *incr_val; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateIncrement() { AKANTU_DEBUG_IN(); - AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." - << " Are you working with Finite or Ineslactic deformations?"); + AKANTU_DEBUG_ASSERT( + previous_displacement, + "The previous displacement has to be initialized." + << " Are you working with Finite or Ineslactic deformations?"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); Real * prev_disp_val = previous_displacement->storage(); - for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++prev_disp_val) + for (UInt j = 0; j < nb_degree_of_freedom; + ++j, ++disp_val, ++incr_val, ++prev_disp_val) *incr_val = (*disp_val - *prev_disp_val); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updatePreviousDisplacement() { AKANTU_DEBUG_IN(); - AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." - << " Are you working with Finite or Ineslactic deformations?"); + AKANTU_DEBUG_ASSERT( + previous_displacement, + "The previous displacement has to be initialized." + << " Are you working with Finite or Ineslactic deformations?"); previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); - AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." - << " Did you call initParallel?"); + AKANTU_DEBUG_ASSERT(synch_registry, + "Synchronizer registry was not initialized." + << " Did you call initParallel?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); - AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." - << " Did you call initPBC?"); + AKANTU_DEBUG_ASSERT(synch_registry, + "Synchronizer registry was not initialized." + << " Did you call initPBC?"); synch_registry->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); - if(!increment) { + if (!increment) { UInt nb_nodes = mesh.getNbNodes(); - std::stringstream sstr_inc; sstr_inc << id << ":increment"; + std::stringstream sstr_inc; + sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0.)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; - Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_regular); - Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type, _ek_regular); - for(; it != end; ++it) { + Mesh::type_iterator it = + mesh.firstType(spatial_dimension, ghost_type, _ek_regular); + Mesh::type_iterator end = + mesh.lastType(spatial_dimension, ghost_type, _ek_regular); + for (; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); - UInt nb_element = mesh.getNbElement(*it); + UInt nb_element = mesh.getNbElement(*it); - Array::const_scalar_iterator mat_indexes = material_index(*it, ghost_type).begin(); - Array::const_scalar_iterator mat_loc_num = material_local_numbering(*it, ghost_type).begin(); + Array::const_scalar_iterator mat_indexes = + material_index(*it, ghost_type).begin(); + Array::const_scalar_iterator mat_loc_num = + material_local_numbering(*it, ghost_type).begin(); - Array X(0, nb_nodes_per_element*spatial_dimension); - FEEngine::extractNodalToElementField(mesh, *current_position, - X, *it, _not_ghost); + Array X(0, nb_nodes_per_element * spatial_dimension); + FEEngine::extractNodalToElementField(mesh, *current_position, X, *it, + _not_ghost); Array::matrix_iterator X_el = - X.begin(spatial_dimension, nb_nodes_per_element); + X.begin(spatial_dimension, nb_nodes_per_element); - for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { + for (UInt el = 0; el < nb_element; + ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; - Real el_h = getFEEngine().getElementInradius(*X_el, *it); - Real el_c = mat_val[*mat_indexes]->getCelerity(elem); + Real el_h = getFEEngine().getElementInradius(*X_el, *it); + Real el_c = mat_val[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); if (!mass && !mass_matrix) AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); - Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); - Real * vel_val = velocity->storage(); + Real * vel_val = velocity->storage(); Real * mass_val = mass->storage(); for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) - mv2 += *vel_val * *vel_val * *mass_val; + mv2 += *vel_val * *vel_val * *mass_val; vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ -Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { +Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, + UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, - spatial_dimension, - type, _not_ghost, - filter_element); + spatial_dimension, type, + _not_ghost, filter_element); - Array::vector_iterator vit = vel_on_quad.begin(spatial_dimension); - Array::vector_iterator vend = vel_on_quad.end(spatial_dimension); + Array::vector_iterator vit = vel_on_quad.begin(spatial_dimension); + Array::vector_iterator vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); - return .5*getFEEngine().integrate(rho_v2, type, index); + return .5 * getFEEngine().integrate(rho_v2, type, index); } - /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * incr_or_velo = NULL; - if(this->method == _static){ + if (this->method == _static) { incr_or_velo = this->increment->storage(); - } - else - incr_or_velo = this->velocity->storage(); + } else + incr_or_velo = this->velocity->storage(); Real * forc = this->force->storage(); Real * resi = this->residual->storage(); bool * boun = this->blocked_dofs->storage(); Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = this->mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < this->spatial_dimension; ++i) { if (count_node) { - if(*boun) - work -= *resi * *incr_or_velo; - else - work += *forc * *incr_or_velo; + if (*boun) + work -= *resi * *incr_or_velo; + else + work += *forc * *incr_or_velo; } ++incr_or_velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); - if(this->method != _static) + if (this->method != _static) work *= this->time_step; AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); - } else if (energy_id == "external work"){ + } else if (energy_id == "external work") { return getExternalWork(); } Real energy = 0.; std::vector::iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getEnergy(energy_id); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, - const ElementType & type, - UInt index) { + const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } std::vector::iterator mat_it; - UInt mat_index = this->material_index(type, _not_ghost)(index); + UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); - Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); + Real energy = + this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, - const NewElementsEvent & event) { + const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); - for(UInt g = _not_ghost; g <= _ghost; ++g) { - GhostType gt = (GhostType) g; - Mesh::type_iterator it = this->mesh.firstType(spatial_dimension, gt, _ek_not_defined); - Mesh::type_iterator end = this->mesh.lastType(spatial_dimension, gt, _ek_not_defined); - for(; it != end; ++it) { + for (UInt g = _not_ghost; g <= _ghost; ++g) { + GhostType gt = (GhostType)g; + Mesh::type_iterator it = + this->mesh.firstType(spatial_dimension, gt, _ek_not_defined); + Mesh::type_iterator end = + this->mesh.lastType(spatial_dimension, gt, _ek_not_defined); + for (; it != end; ++it) { UInt nb_element = this->mesh.getNbElement(*it, gt); - if(!material_index.exists(*it, gt)) { - this->material_index .alloc(nb_element, 1, *it, gt); - this->material_local_numbering.alloc(nb_element, 1, *it, gt); + if (!material_index.exists(*it, gt)) { + this->material_index.alloc(nb_element, 1, *it, gt); + this->material_local_numbering.alloc(nb_element, 1, *it, gt); } else { - this->material_index (*it, gt).resize(nb_element); - this->material_local_numbering(*it, gt).resize(nb_element); + this->material_index(*it, gt).resize(nb_element); + this->material_local_numbering(*it, gt).resize(nb_element); } } } - Array::const_iterator it = element_list.begin(); + Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); ElementTypeMapArray filter("new_element_filter", this->getID()); for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; - if(!filter.exists(elem.type, elem.ghost_type)) + if (!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); std::vector::iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsAdded(element_list, event); } - if(method == _explicit_lumped_mass) this->assembleMassLumped(); + if (method == _explicit_lumped_mass) + this->assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array & element_list, - const ElementTypeMapArray & new_numbering, - const RemovedElementsEvent & event) { +void SolidMechanicsModel::onElementsRemoved( + __attribute__((unused)) const Array & element_list, + const ElementTypeMapArray & new_numbering, + const RemovedElementsEvent & event) { this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); std::vector::iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering, event); } - } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, - __attribute__((unused)) const NewNodesEvent & event) { + __attribute__((unused)) + const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); - if(displacement) displacement->resize(nb_nodes); - if(mass ) mass ->resize(nb_nodes); - if(velocity ) velocity ->resize(nb_nodes); - if(acceleration) acceleration->resize(nb_nodes); - if(force ) force ->resize(nb_nodes); - if(residual ) residual ->resize(nb_nodes); - if(blocked_dofs) blocked_dofs->resize(nb_nodes); - - if(previous_displacement) previous_displacement->resize(nb_nodes); - if(increment_acceleration) increment_acceleration->resize(nb_nodes); - if(increment) increment->resize(nb_nodes); - - if(current_position) current_position->resize(nb_nodes); + if (displacement) + displacement->resize(nb_nodes); + if (mass) + mass->resize(nb_nodes); + if (velocity) + velocity->resize(nb_nodes); + if (acceleration) + acceleration->resize(nb_nodes); + if (force) + force->resize(nb_nodes); + if (residual) + residual->resize(nb_nodes); + if (blocked_dofs) + blocked_dofs->resize(nb_nodes); + + if (previous_displacement) + previous_displacement->resize(nb_nodes); + if (increment_acceleration) + increment_acceleration->resize(nb_nodes); + if (increment) + increment->resize(nb_nodes); + + if (current_position) + current_position->resize(nb_nodes); std::vector::iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); if (method != _explicit_lumped_mass) { this->initSolver(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, - const Array & new_numbering, - __attribute__((unused)) const RemovedNodesEvent & event) { - if(displacement) mesh.removeNodesFromArray(*displacement, new_numbering); - if(mass ) mesh.removeNodesFromArray(*mass , new_numbering); - if(velocity ) mesh.removeNodesFromArray(*velocity , new_numbering); - if(acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); - if(force ) mesh.removeNodesFromArray(*force , new_numbering); - if(residual ) mesh.removeNodesFromArray(*residual , new_numbering); - if(blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); - - if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering); - if(increment) mesh.removeNodesFromArray(*increment , new_numbering); +void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) + const Array & element_list, + const Array & new_numbering, + __attribute__((unused)) + const RemovedNodesEvent & event) { + if (displacement) + mesh.removeNodesFromArray(*displacement, new_numbering); + if (mass) + mesh.removeNodesFromArray(*mass, new_numbering); + if (velocity) + mesh.removeNodesFromArray(*velocity, new_numbering); + if (acceleration) + mesh.removeNodesFromArray(*acceleration, new_numbering); + if (force) + mesh.removeNodesFromArray(*force, new_numbering); + if (residual) + mesh.removeNodesFromArray(*residual, new_numbering); + if (blocked_dofs) + mesh.removeNodesFromArray(*blocked_dofs, new_numbering); + + if (increment_acceleration) + mesh.removeNodesFromArray(*increment_acceleration, new_numbering); + if (increment) + mesh.removeNodesFromArray(*increment, new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); if (method != _explicit_lumped_mass) { this->initSolver(); } - } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::isInternal(const std::string & field_name, const ElementKind & element_kind) { bool is_internal = false; /// check if at least one material contains field_id as an internal for (UInt m = 0; m < materials.size() && !is_internal; ++m) { is_internal |= materials[m]->isInternal(field_name, element_kind); } return is_internal; } /* -------------------------------------------------------------------------- */ ElementTypeMap SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name, const ElementKind & element_kind) { if (!(this->isInternal(field_name, element_kind))) AKANTU_EXCEPTION("unknown internal " << field_name); for (UInt m = 0; m < materials.size(); ++m) { if (materials[m]->isInternal(field_name, element_kind)) - return materials[m]->getInternalDataPerElem(field_name, element_kind); + return materials[m]->getInternalDataPerElem(field_name, + element_kind); } return ElementTypeMap(); } /* -------------------------------------------------------------------------- */ ElementTypeMapArray & SolidMechanicsModel::flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type) { std::pair key(field_name, kind); if (this->registered_internals.count(key) == 0) { this->registered_internals[key] = new ElementTypeMapArray(field_name, this->id); } ElementTypeMapArray * internal_flat = this->registered_internals[key]; typedef ElementTypeMapArray::type_iterator iterator; iterator tit = internal_flat->firstType(spatial_dimension, ghost_type, kind); iterator end = internal_flat->lastType(spatial_dimension, ghost_type, kind); for (; tit != end; ++tit) { ElementType type = *tit; (*internal_flat)(type, ghost_type).clear(); } for (UInt m = 0; m < materials.size(); ++m) { if (materials[m]->isInternal(field_name, kind)) materials[m]->flattenInternal(field_name, *internal_flat, ghost_type, kind); } return *internal_flat; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::flattenAllRegisteredInternals( const ElementKind & kind) { typedef std::map, ElementTypeMapArray *>::iterator iterator; iterator it = this->registered_internals.begin(); iterator end = this->registered_internals.end(); while (it != end) { if (kind == it->first.second) this->flattenInternal(it->first.first, kind); ++it; } } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModel::onDump(){ +void SolidMechanicsModel::onDump() { this->flattenAllRegisteredInternals(_ek_regular); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER -dumper::Field * SolidMechanicsModel -::createElementalField(const std::string & field_name, - const std::string & group_name, - bool padding_flag, - const UInt & spatial_dimension, - const ElementKind & kind) { +dumper::Field * SolidMechanicsModel::createElementalField( + const std::string & field_name, const std::string & group_name, + bool padding_flag, const UInt & spatial_dimension, + const ElementKind & kind) { dumper::Field * field = NULL; - if(field_name == "partitions") - field = mesh.createElementalField(mesh.getConnectivities(),group_name,spatial_dimension,kind); - else if(field_name == "material_index") - field = mesh.createElementalField(material_index,group_name,spatial_dimension,kind); + if (field_name == "partitions") + field = mesh.createElementalField( + mesh.getConnectivities(), group_name, spatial_dimension, kind); + else if (field_name == "material_index") + field = mesh.createElementalField( + material_index, group_name, spatial_dimension, kind); else { // this copy of field_name is used to compute derivated data such as // strain and von mises stress that are based on grad_u and stress std::string field_name_copy(field_name); - if (field_name == "strain" - || field_name == "Green strain" - || field_name == "principal strain" - || field_name == "principal Green strain") + if (field_name == "strain" || field_name == "Green strain" || + field_name == "principal strain" || + field_name == "principal Green strain") field_name_copy = "grad_u"; else if (field_name == "Von Mises stress") field_name_copy = "stress"; - bool is_internal = this->isInternal(field_name_copy,kind); + bool is_internal = this->isInternal(field_name_copy, kind); if (is_internal) { - ElementTypeMap nb_data_per_elem = this->getInternalDataPerElem(field_name_copy, kind); - ElementTypeMapArray & internal_flat = this->flattenInternal(field_name_copy,kind); - field = mesh.createElementalField(internal_flat, - group_name, - spatial_dimension,kind,nb_data_per_elem); - if (field_name == "strain"){ - dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); - field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); + ElementTypeMap nb_data_per_elem = + this->getInternalDataPerElem(field_name_copy, kind); + ElementTypeMapArray & internal_flat = + this->flattenInternal(field_name_copy, kind); + field = mesh.createElementalField( + internal_flat, group_name, spatial_dimension, kind, nb_data_per_elem); + if (field_name == "strain") { + dumper::ComputeStrain * foo = + new dumper::ComputeStrain(*this); + field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "Von Mises stress") { - dumper::ComputeVonMisesStress * foo = new dumper::ComputeVonMisesStress(*this); - field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); + dumper::ComputeVonMisesStress * foo = + new dumper::ComputeVonMisesStress(*this); + field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "Green strain") { - dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); - field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); + dumper::ComputeStrain * foo = + new dumper::ComputeStrain(*this); + field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "principal strain") { - dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); - field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); + dumper::ComputePrincipalStrain * foo = + new dumper::ComputePrincipalStrain(*this); + field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } else if (field_name == "principal Green strain") { - dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); - field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); + dumper::ComputePrincipalStrain * foo = + new dumper::ComputePrincipalStrain(*this); + field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } - //treat the paddings - if (padding_flag){ - if (field_name == "stress"){ - if (spatial_dimension == 2) { - dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this); - field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); - } - } else if (field_name == "strain" || field_name == "Green strain"){ - if (spatial_dimension == 2) { - dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this); - field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); - } - } + // treat the paddings + if (padding_flag) { + if (field_name == "stress") { + if (spatial_dimension == 2) { + dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this); + field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); + } + } else if (field_name == "strain" || field_name == "Green strain") { + if (spatial_dimension == 2) { + dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this); + field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); + } + } } // homogenize the field dumper::ComputeFunctorInterface * foo = - dumper::HomogenizerProxy::createHomogenizer(*field); - - field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); + dumper::HomogenizerProxy::createHomogenizer(*field); + field = dumper::FieldComputeProxy::createFieldCompute(field, *foo); } } return field; } /* -------------------------------------------------------------------------- */ -dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag) { +dumper::Field * +SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, + const std::string & group_name, + bool padding_flag) { - std::map* > real_nodal_fields; - real_nodal_fields["displacement" ] = displacement; - real_nodal_fields["mass" ] = mass; - real_nodal_fields["velocity" ] = velocity; - real_nodal_fields["acceleration" ] = acceleration; - real_nodal_fields["force" ] = force; - real_nodal_fields["residual" ] = residual; - real_nodal_fields["increment" ] = increment; + std::map *> real_nodal_fields; + real_nodal_fields["displacement"] = displacement; + real_nodal_fields["mass"] = mass; + real_nodal_fields["velocity"] = velocity; + real_nodal_fields["acceleration"] = acceleration; + real_nodal_fields["force"] = force; + real_nodal_fields["residual"] = residual; + real_nodal_fields["increment"] = increment; dumper::Field * field = NULL; if (padding_flag) field = mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); else field = mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ -dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, - const std::string & group_name, -bool padding_flag) { +dumper::Field * +SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, + const std::string & group_name, + bool padding_flag) { - std::map* > uint_nodal_fields; - uint_nodal_fields["blocked_dofs" ] = blocked_dofs; + std::map *> uint_nodal_fields; + uint_nodal_fields["blocked_dofs"] = blocked_dofs; dumper::Field * field = NULL; - field = mesh.createNodalField(uint_nodal_fields[field_name],group_name); + field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; - } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ -dumper::Field * SolidMechanicsModel -::createElementalField(const std::string & field_name, - const std::string & group_name, - bool padding_flag, - const UInt & spatial_dimension, - const ElementKind & kind){ +dumper::Field * SolidMechanicsModel::createElementalField( + const std::string & field_name, const std::string & group_name, + bool padding_flag, const UInt & spatial_dimension, + const ElementKind & kind) { return NULL; } /* -------------------------------------------------------------------------- */ -dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag) { +dumper::Field * +SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, + const std::string & group_name, + bool padding_flag) { return NULL; } /* -------------------------------------------------------------------------- */ -dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, - const std::string & group_name, -bool padding_flag) { +dumper::Field * +SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, + const std::string & group_name, + bool padding_flag) { return NULL; } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ -void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { +void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, + UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump() { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeCauchyStresses() { AKANTU_DEBUG_IN(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; - if(mat.isFiniteDeformation()) + if (mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::saveStressAndStrainBeforeDamage() { - EventManager::sendEvent(SolidMechanicsModelEvent::BeginningOfDamageIterationEvent()); + EventManager::sendEvent( + SolidMechanicsModelEvent::BeginningOfDamageIterationEvent()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateEnergiesAfterDamage() { EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; - for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); + for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) + ; stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; - stream << space << " + spatial dimension : " << spatial_dimension << std::endl; + stream << space << " + spatial dimension : " << spatial_dimension + << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); - mass ->printself(stream, indent + 2); - velocity ->printself(stream, indent + 2); + mass->printself(stream, indent + 2); + velocity->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); - force ->printself(stream, indent + 2); - residual ->printself(stream, indent + 2); + force->printself(stream, indent + 2); + residual->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; std::vector::const_iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { const Material & mat = *(*mat_it); mat.printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 4595f33ae..aeeae6b85 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,775 +1,780 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Dec 18 2015 * * @brief Model of Solid Mechanics * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ - /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" #include "model.hh" #include "data_accessor.hh" #include "mesh.hh" #include "dumpable.hh" #include "boundary_condition.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "integration_scheme_2nd_order.hh" #include "solver.hh" #include "material_selector.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ namespace akantu { - class Material; - class IntegrationScheme2ndOrder; - class SparseMatrix; - class DumperIOHelper; - class NonLocalManager; +class Material; +class IntegrationScheme2ndOrder; +class SparseMatrix; +class DumperIOHelper; +class NonLocalManager; } /* -------------------------------------------------------------------------- */ - __BEGIN_AKANTU__ struct SolidMechanicsModelOptions : public ModelOptions { - SolidMechanicsModelOptions(AnalysisMethod analysis_method = _explicit_lumped_mass, - bool no_init_materials = false) : - analysis_method(analysis_method), - no_init_materials(no_init_materials) { } + SolidMechanicsModelOptions( + AnalysisMethod analysis_method = _explicit_lumped_mass, + bool no_init_materials = false) + : analysis_method(analysis_method), no_init_materials(no_init_materials) { + } AnalysisMethod analysis_method; bool no_init_materials; }; extern const SolidMechanicsModelOptions default_solid_mechanics_model_options; - -class SolidMechanicsModel : public Model, - public DataAccessor, - public MeshEventHandler, - public BoundaryCondition, - public EventHandlerManager { +class SolidMechanicsModel + : public Model, + public DataAccessor, + public MeshEventHandler, + public BoundaryCondition, + public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: - AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); - AKANTU_GET_MACRO(MaterialList, material, const Array &); + AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); + AKANTU_GET_MACRO(MaterialList, material, const Array &); + protected: - Array material; + Array material; }; - typedef FEEngineTemplate MyFEEngineType; + typedef FEEngineTemplate MyFEEngineType; protected: - typedef EventHandlerManager EventManager; + typedef EventHandlerManager EventManager; public: - SolidMechanicsModel(Mesh & mesh, - UInt spatial_dimension = _all_dimensions, - const ID & id = "solid_mechanics_model", - const MemoryID & memory_id = 0); + SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, + const ID & id = "solid_mechanics_model", + const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model - virtual void initFull(const ModelOptions & options = default_solid_mechanics_model_options); + virtual void initFull( + const ModelOptions & options = default_solid_mechanics_model_options); /// initialize the fem object needed for boundary conditions void initFEEngineBoundary(); /// register the tags associated with the parallel synchronizer - virtual void initParallel(MeshPartition *partition, - DataAccessor *data_accessor = NULL); + virtual void initParallel(MeshPartition * partition, + DataAccessor * data_accessor = NULL); /// allocate all vectors virtual void initArrays(); /// allocate all vectors void initArraysPreviousDisplacment(); /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model virtual void initModel(); /// init PBC synchronizer void initPBC(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// re-initialize model to set fields to 0 void reInitialize(); /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC // void changeEquationNumberforPBC(std::map & pbc_pair); /// synchronize Residual for output void synchronizeResidual(); protected: /// register PBC synchronizer void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the stuff for the explicit scheme void initExplicit(AnalysisMethod analysis_method = _explicit_lumped_mass); bool isExplicit() { - return method == _explicit_lumped_mass || method == _explicit_consistent_mass; + return method == _explicit_lumped_mass || + method == _explicit_consistent_mass; } /// initialize the array needed by updateResidual (residual, current_position) void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); /// assemble the residual for the explicit scheme virtual void updateResidual(bool need_initialize = true); /** * \brief compute the acceleration from the residual * this function is the explicit equivalent to solveDynamic in implicit * In the case of lumped mass just divide the residual by the mass * In the case of not lumped mass call solveDynamic<_acceleration_corrector> */ void updateAcceleration(); - ///Update the increment of displacement + /// Update the increment of displacement void updateIncrement(); - ///Copy the actuel displacement into previous displacement + /// Copy the actuel displacement into previous displacement void updatePreviousDisplacement(); - ///Save stress and strain through EventManager + /// Save stress and strain through EventManager void saveStressAndStrainBeforeDamage(); - ///Update energies through EventManager + /// Update energies through EventManager void updateEnergiesAfterDamage(); /// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix - void solveLumped(Array & x, - const Array & A, - const Array & b, - const Array & blocked_dofs, - Real alpha); + void solveLumped(Array & x, const Array & A, + const Array & b, const Array & blocked_dofs, + Real alpha); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); public: void solveStep(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the solver and the jacobian_matrix (called by initImplicit) void initSolver(SolverOptions & options = _solver_no_options); /// initialize the stuff for the implicit solver - void initImplicit(bool dynamic = false, - SolverOptions & solver_options = _solver_no_options); + void initImplicit(bool dynamic = false, + SolverOptions & solver_options = _solver_no_options); /// solve Ma = f to get the initial acceleration void initialAcceleration(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); public: /** * solve a step (predictor + convergence loop + corrector) using the * the given convergence method (see akantu::SolveConvergenceMethod) * and the given convergence criteria (see * akantu::SolveConvergenceCriteria) **/ template bool solveStep(Real tolerance, UInt max_iteration = 100); template - bool solveStep(Real tolerance, - Real & error, - UInt max_iteration = 100, - bool do_not_factorize = false); + bool solveStep(Real tolerance, Real & error, UInt max_iteration = 100, + bool do_not_factorize = false); public: /** * solve Ku = f using the the given convergence method (see * akantu::SolveConvergenceMethod) and the given convergence * criteria (see akantu::SolveConvergenceCriteria) **/ template bool solveStatic(Real tolerance, UInt max_iteration, - bool do_not_factorize = false); + bool do_not_factorize = false); /// test if the system is converged template bool testConvergence(Real tolerance, Real & error); /// test the convergence (norm of increment) bool testConvergenceIncrement(Real tolerance) __attribute__((deprecated)); - bool testConvergenceIncrement(Real tolerance, Real & error) __attribute__((deprecated)); + bool testConvergenceIncrement(Real tolerance, Real & error) + __attribute__((deprecated)); /// test the convergence (norm of residual) bool testConvergenceResidual(Real tolerance) __attribute__((deprecated)); - bool testConvergenceResidual(Real tolerance, Real & error) __attribute__((deprecated)); + bool testConvergenceResidual(Real tolerance, Real & error) + __attribute__((deprecated)); /// create and return the velocity damping matrix SparseMatrix & initVelocityDampingMatrix(); /// implicit time integration predictor void implicitPred(); /// implicit time integration corrector void implicitCorr(); /// compute the Cauchy stress on user demand. void computeCauchyStresses(); /// compute A and solve @f[ A\delta u = f_ext - f_int @f] template - void solve(Array &increment, Real block_val = 1., + void solve(Array & increment, Real block_val = 1., bool need_factorize = true, bool has_profile_changed = false); protected: /// finish the computation of residual to solve in increment void updateResidualInternal(); /// compute the support reaction and store it in force void updateSupportReaction(); private: /// re-initialize the J matrix (to use if the profile of K changed) void initJacobianMatrix(); /* ------------------------------------------------------------------------ */ /* Explicit/Implicit */ /* ------------------------------------------------------------------------ */ public: /// Update the stresses for the computation of the residual of the Stiffness /// matrix in the case of finite deformation void computeStresses(); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: - /// registers all the custom materials of a given type present in the input file - template - void registerNewCustomMaterials(const ID & mat_type); + /// registers all the custom materials of a given type present in the input + /// file + template void registerNewCustomMaterials(const ID & mat_type); /// register an empty material of a given type template Material & registerNewEmptyMaterial(const std::string & name); // /// Use a UIntData in the mesh to specify the material to use per element // void setMaterialIDsFromIntData(const std::string & data_name); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material - virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); - + virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, + const ID & material_name, + const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database template Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials - void assignMaterialToElements(const ElementTypeMapArray * filter = NULL); + void + assignMaterialToElements(const ElementTypeMapArray * filter = NULL); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); - protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho - void computeRho(Array & rho, - ElementType type, - GhostType ghost_type); + void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); - /// compute the external work (for impose displacement, the velocity should be given too) + /// compute the external work (for impose displacement, the velocity should be + /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: + inline virtual UInt getNbDataForElements(const Array & elements, + SynchronizationTag tag) const; - inline virtual UInt getNbDataForElements(const Array & elements, - SynchronizationTag tag) const; + inline virtual void packElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag) const; - inline virtual void packElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag) const; - - inline virtual void unpackElementData(CommunicationBuffer & buffer, - const Array & elements, - SynchronizationTag tag); + inline virtual void unpackElementData(CommunicationBuffer & buffer, + const Array & elements, + SynchronizationTag tag); inline virtual UInt getNbDataToPack(SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; - inline virtual void packData(CommunicationBuffer & buffer, - const UInt index, - SynchronizationTag tag) const; + inline virtual void packData(CommunicationBuffer & buffer, const UInt index, + SynchronizationTag tag) const; - inline virtual void unpackData(CommunicationBuffer & buffer, - const UInt index, - SynchronizationTag tag); + inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, + SynchronizationTag tag); protected: - inline void splitElementByMaterial(const Array & elements, - Array * elements_per_mat) const; + inline void splitElementByMaterial(const Array & elements, + Array * elements_per_mat) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: - virtual void onNodesAdded(const Array & nodes_list, - const NewNodesEvent & event); - virtual void onNodesRemoved(const Array & element_list, - const Array & new_numbering, - const RemovedNodesEvent & event); - virtual void onElementsAdded(const Array & nodes_list, - const NewElementsEvent & event); - virtual void onElementsRemoved(const Array & element_list, - const ElementTypeMapArray & new_numbering, - const RemovedElementsEvent & event); - - virtual void onElementsChanged(__attribute__((unused)) const Array & old_elements_list, - __attribute__((unused)) const Array & new_elements_list, - __attribute__((unused)) const ElementTypeMapArray & new_numbering, - __attribute__((unused)) const ChangedElementsEvent & event) {}; - + virtual void onNodesAdded(const Array & nodes_list, + const NewNodesEvent & event); + virtual void onNodesRemoved(const Array & element_list, + const Array & new_numbering, + const RemovedNodesEvent & event); + virtual void onElementsAdded(const Array & nodes_list, + const NewElementsEvent & event); + virtual void + onElementsRemoved(const Array & element_list, + const ElementTypeMapArray & new_numbering, + const RemovedElementsEvent & event); + + virtual void onElementsChanged( + __attribute__((unused)) const Array & old_elements_list, + __attribute__((unused)) const Array & new_elements_list, + __attribute__((unused)) const ElementTypeMapArray & new_numbering, + __attribute__((unused)) const ChangedElementsEvent & event){}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: - - virtual void onDump(); //! decide wether a field is a material internal or not - bool isInternal(const std::string & field_name, const ElementKind & element_kind); + bool isInternal(const std::string & field_name, + const ElementKind & element_kind); #ifndef SWIG //! give the amount of data per element - virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, - const ElementKind & kind); + virtual ElementTypeMap + getInternalDataPerElem(const std::string & field_name, + const ElementKind & kind); //! flatten a given material internal field - ElementTypeMapArray & flattenInternal(const std::string & field_name, - const ElementKind & kind, - const GhostType ghost_type = _not_ghost); + ElementTypeMapArray & + flattenInternal(const std::string & field_name, const ElementKind & kind, + const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); #endif virtual dumper::Field * createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag); + const std::string & group_name, + bool padding_flag); virtual dumper::Field * createNodalFieldBool(const std::string & field_name, - const std::string & group_name, - bool padding_flag); - + const std::string & group_name, + bool padding_flag); virtual dumper::Field * createElementalField(const std::string & field_name, - const std::string & group_name, - bool padding_flag, - const UInt & spatial_dimension, - const ElementKind & kind); + const std::string & group_name, + bool padding_flag, + const UInt & spatial_dimension, + const ElementKind & kind); virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); virtual void dump(); virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step void setTimeStep(Real time_step); /// return the of iterations done in the last solveStep AKANTU_GET_MACRO(NumberIter, n_iter, UInt); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector - AKANTU_GET_MACRO(Displacement, *displacement, Array &); + AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector - AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); + AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition - AKANTU_GET_MACRO(CurrentPosition, *current_position, const Array &); + AKANTU_GET_MACRO(CurrentPosition, *current_position, const Array &); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before - AKANTU_GET_MACRO(Increment, *increment, Array &); + AKANTU_GET_MACRO(Increment, *increment, Array &); /// get the lumped SolidMechanicsModel::mass vector - AKANTU_GET_MACRO(Mass, *mass, Array &); + AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector - AKANTU_GET_MACRO(Velocity, *velocity, Array &); + AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration - AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); + AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::force vector (boundary forces) - AKANTU_GET_MACRO(Force, *force, Array &); + AKANTU_GET_MACRO(Force, *force, Array &); /// get the SolidMechanicsModel::residual vector, computed by /// SolidMechanicsModel::updateResidual - AKANTU_GET_MACRO(Residual, *residual, Array &); + AKANTU_GET_MACRO(Residual, *residual, Array &); /// get the SolidMechanicsModel::blocked_dofs vector - AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); + AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the SolidMechnicsModel::incrementAcceleration vector - AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, Array &); + AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, + Array &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials - inline UInt getNbMaterials() const { - return materials.size(); - } + inline UInt getNbMaterials() const { return materials.size(); } inline void setMaterialSelector(MaterialSelector & selector); /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy - Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); + Real getEnergy(const std::string & energy_id, const ElementType & type, + UInt index); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ void setIncrementFlagOn(); /// get the stiffness matrix AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); /// get the global jacobian matrix of the system - AKANTU_GET_MACRO(GlobalJacobianMatrix, *jacobian_matrix, const SparseMatrix &); + AKANTU_GET_MACRO(GlobalJacobianMatrix, *jacobian_matrix, + const SparseMatrix &); /// get the mass matrix AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); /// get the velocity damping matrix - AKANTU_GET_MACRO(VelocityDampingMatrix, *velocity_damping_matrix, SparseMatrix &); + AKANTU_GET_MACRO(VelocityDampingMatrix, *velocity_damping_matrix, + SparseMatrix &); /// get the FEEngine object to integrate or interpolate on the boundary inline FEEngine & getFEEngineBoundary(const ID & name = ""); /// get integrator AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder &); /// get access to the internal solver AKANTU_GET_MACRO(Solver, *solver, Solver &); /// get synchronizer - AKANTU_GET_MACRO(Synchronizer, *synch_parallel, const DistributedSynchronizer &); + AKANTU_GET_MACRO(Synchronizer, *synch_parallel, + const DistributedSynchronizer &); - AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); + AKANTU_GET_MACRO(MaterialByElement, material_index, + const ElementTypeMapArray &); - /// vectors containing local material element index for each global element index - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); + /// vectors containing local material element index for each global element + /// index + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, + UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); - AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); + AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, + material_local_numbering, UInt); + AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, + material_local_numbering, UInt); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /// get the non-local manager AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); - template - friend struct ContactData; + template friend struct ContactData; template friend class ContactResolution; protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// number of iterations UInt n_iter; /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array - Array *displacement; + Array * displacement; /// displacements array at the previous time step (used in finite deformation) - Array *previous_displacement; + Array * previous_displacement; /// lumped mass array - Array *mass; + Array * mass; /// velocities array - Array *velocity; + Array * velocity; /// accelerations array - Array *acceleration; + Array * acceleration; /// accelerations array - Array *increment_acceleration; + Array * increment_acceleration; /// forces array - Array *force; + Array * force; /// residuals array - Array *residual; + Array * residual; /// array specifing if a degree of freedom is blocked or not - Array *blocked_dofs; + Array * blocked_dofs; /// array of current position used during update residual - Array *current_position; + Array * current_position; /// mass matrix - SparseMatrix *mass_matrix; + SparseMatrix * mass_matrix; /// velocity damping matrix - SparseMatrix *velocity_damping_matrix; + SparseMatrix * velocity_damping_matrix; /// stiffness matrix - SparseMatrix *stiffness_matrix; + SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta /// t^2}, d = \frac{\gamma}{\beta \Delta t} @f] - SparseMatrix *jacobian_matrix; + SparseMatrix * jacobian_matrix; /// Arrays containing the material index for each element ElementTypeMapArray material_index; - /// Arrays containing the position in the element filter of the material (material's local numbering) + /// Arrays containing the position in the element filter of the material + /// (material's local numbering) ElementTypeMapArray material_local_numbering; #ifdef SWIGPYTHON public: #endif /// list of used materials - std::vector materials; + std::vector materials; /// mapping between material name and material internal id - std::map materials_names_to_id; + std::map materials_names_to_id; #ifdef SWIGPYTHON protected: #endif /// class defining of to choose a material - MaterialSelector *material_selector; + MaterialSelector * material_selector; /// define if it is the default selector or not bool is_default_material_selector; /// integration scheme of second order used - IntegrationScheme2ndOrder *integrator; + IntegrationScheme2ndOrder * integrator; /// increment of displacement - Array *increment; + Array * increment; /// flag defining if the increment must be computed or not bool increment_flag; /// solver for implicit - Solver *solver; + Solver * solver; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// internal synchronizer for parallel computations - DistributedSynchronizer *synch_parallel; + DistributedSynchronizer * synch_parallel; /// tells if the material are instantiated bool are_materials_instantiated; + typedef std::map, + ElementTypeMapArray *> flatten_internal_map; + /// map a registered internals to be flattened for dump purposes - std::map,ElementTypeMapArray *> registered_internals; + flatten_internal_map registered_internals; /// pointer to non-local manager: For non-local continuum damage computations - NonLocalManager * non_local_manager; + NonLocalManager * non_local_manager; /// pointer to the pbc synchronizer PBCSynchronizer * pbc_synch; - }; - /* -------------------------------------------------------------------------- */ namespace BC { - namespace Neumann { - typedef FromHigherDim FromStress; - typedef FromSameDim FromTraction; - } +namespace Neumann { + typedef FromHigherDim FromStress; + typedef FromSameDim FromTraction; +} } __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "parser.hh" #include "material.hh" __BEGIN_AKANTU__ #include "solid_mechanics_model_tmpl.hh" -#if defined (AKANTU_INCLUDE_INLINE_IMPL) -# include "solid_mechanics_model_inline_impl.cc" +#if defined(AKANTU_INCLUDE_INLINE_IMPL) +#include "solid_mechanics_model_inline_impl.cc" #endif /// standard output stream operator -inline std::ostream & operator << (std::ostream & stream, const SolidMechanicsModel &_this) { +inline std::ostream & operator<<(std::ostream & stream, + const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ - #include "material_selector_tmpl.hh" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/solver/solver.cc b/src/solver/solver.cc index 5fa0afd38..9c911a805 100644 --- a/src/solver/solver.cc +++ b/src/solver/solver.cc @@ -1,90 +1,86 @@ /** * @file solver.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Wed Oct 07 2015 * * @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 "solver.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ SolverOptions _solver_no_options(true); /* -------------------------------------------------------------------------- */ -Solver::Solver(SparseMatrix & matrix, - const ID & id, - const MemoryID & memory_id) : - Memory(id, memory_id), StaticSolverEventHandler(), - matrix(&matrix), - is_matrix_allocated(false), - mesh(NULL), - communicator(StaticCommunicator::getStaticCommunicator()), - solution(NULL), - synch_registry(NULL) { +Solver::Solver(SparseMatrix & matrix, const ID & id, const MemoryID & memory_id) + : Memory(id, memory_id), StaticSolverEventHandler(), matrix(&matrix), + is_matrix_allocated(false), mesh(NULL), + communicator(StaticCommunicator::getStaticCommunicator()), solution(NULL), + synch_registry(NULL) { AKANTU_DEBUG_IN(); StaticSolver::getStaticSolver().registerEventHandler(*this); - //createSynchronizerRegistry(); + // createSynchronizerRegistry(); this->synch_registry = new SynchronizerRegistry(*this); - synch_registry->registerSynchronizer(this->matrix->getDOFSynchronizer(), _gst_solver_solution); + synch_registry->registerSynchronizer(this->matrix->getDOFSynchronizer(), + _gst_solver_solution); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Solver::~Solver() { AKANTU_DEBUG_IN(); this->destroyInternalData(); delete synch_registry; StaticSolver::getStaticSolver().unregisterEventHandler(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Solver::beforeStaticSolverDestroy() { AKANTU_DEBUG_IN(); - try{ + try { this->destroyInternalData(); - } catch(...) {} + } catch (...) { + } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Solver::createSynchronizerRegistry() { - //this->synch_registry = new SynchronizerRegistry(this); + // this->synch_registry = new SynchronizerRegistry(this); } - __END_AKANTU__ diff --git a/src/solver/solver.hh b/src/solver/solver.hh index dc93acadb..cfb485fa3 100644 --- a/src/solver/solver.hh +++ b/src/solver/solver.hh @@ -1,168 +1,161 @@ /** * @file solver.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Dec 08 2015 * * @brief interface for solvers * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_HH__ #define __AKANTU_SOLVER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_array.hh" #include "sparse_matrix.hh" #include "mesh.hh" #include "static_communicator.hh" #include "static_solver.hh" #include "data_accessor.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class SolverOptions { public: SolverOptions(bool no_option = false) // : no_option(no_option) - { } + {} virtual ~SolverOptions() {} private: - //bool no_option; + // bool no_option; }; extern SolverOptions _solver_no_options; class Solver : protected Memory, - public StaticSolverEventHandler, - public DataAccessor { + public StaticSolverEventHandler, + public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - - Solver(SparseMatrix & matrix, - const ID & id = "solver", - const MemoryID & memory_id = 0); + Solver(SparseMatrix & matrix, const ID & id = "solver", + const MemoryID & memory_id = 0); virtual ~Solver(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// initialize the solver virtual void initialize(SolverOptions & options = _solver_no_options) = 0; - virtual void setOperators() {}; - virtual void analysis() {}; + virtual void setOperators(){}; + virtual void analysis(){}; - virtual void factorize() {}; + virtual void factorize(){}; /// solve virtual void solve(Array & solution) = 0; virtual void solve() = 0; virtual void setRHS(Array & rhs) = 0; /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const; protected: - virtual void destroyInternalData() {}; + virtual void destroyInternalData(){}; public: virtual void beforeStaticSolverDestroy(); void createSynchronizerRegistry(); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: - inline virtual UInt getNbDataForDOFs(const Array & dofs, - SynchronizationTag tag) const; + inline virtual UInt getNbDataForDOFs(const Array & dofs, + SynchronizationTag tag) const; inline virtual void packDOFData(CommunicationBuffer & buffer, - const Array & dofs, - SynchronizationTag tag) const; + const Array & dofs, + SynchronizationTag tag) const; inline virtual void unpackDOFData(CommunicationBuffer & buffer, - const Array & dofs, - SynchronizationTag tag); + const Array & dofs, + SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - AKANTU_GET_MACRO(RHS, *rhs, Array &); - /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - /// the matrix SparseMatrix * matrix; /// is sparse matrix allocated bool is_matrix_allocated; /// the right hand side Array * rhs; /// mesh const Mesh * mesh; /// pointer to the communicator StaticCommunicator & communicator; /// the solution obtained from the solve step Array * solution; /// synchronizer registry SynchronizerRegistry * synch_registry; - }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "solver_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_SOLVER_HH__ */ diff --git a/src/solver/solver_mumps.cc b/src/solver/solver_mumps.cc index 9e1055037..6c8e2295b 100644 --- a/src/solver/solver_mumps.cc +++ b/src/solver/solver_mumps.cc @@ -1,390 +1,400 @@ /** * @file solver_mumps.cc * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Fri Oct 16 2015 * * @brief implem of SolverMumps 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" #if defined(AKANTU_USE_MPI) -# include "static_communicator_mpi.hh" -# include "mpi_type_wrapper.hh" +#include "static_communicator_mpi.hh" +#include "mpi_type_wrapper.hh" #endif #include "solver_mumps.hh" #include "dof_synchronizer.hh" - /* -------------------------------------------------------------------------- */ -// static std::ostream & operator <<(std::ostream & stream, const DMUMPS_STRUC_C & _this) { +// 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; // } __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -SolverMumps::SolverMumps(SparseMatrix & matrix, - const ID & id, - const MemoryID & memory_id) : - Solver(matrix, id, memory_id), is_mumps_data_initialized(false), rhs_is_local(true) { +SolverMumps::SolverMumps(SparseMatrix & matrix, const ID & id, + const MemoryID & memory_id) + : Solver(matrix, id, memory_id), is_mumps_data_initialized(false), + rhs_is_local(true) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_MPI parallel_method = SolverMumpsOptions::_fully_distributed; -#else //AKANTU_USE_MPI +#else // AKANTU_USE_MPI parallel_method = SolverMumpsOptions::_not_parallel; -#endif //AKANTU_USE_MPI +#endif // AKANTU_USE_MPI AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolverMumps::~SolverMumps() { AKANTU_DEBUG_IN(); this->destroyInternalData(); - + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::destroyInternalData() { AKANTU_DEBUG_IN(); - if(this->is_mumps_data_initialized) { + if (this->is_mumps_data_initialized) { this->mumps_data.job = _smj_destroy; // destroy dmumps_c(&this->mumps_data); this->is_mumps_data_initialized = false; } AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ void SolverMumps::initMumpsData() { // Default Scaling icntl(8) = 77; // Assembled matrix icntl(5) = 0; /// Default centralized dense second member icntl(20) = 0; icntl(21) = 0; - icntl(28) = 0; //automatic choice for analysis analysis + icntl(28) = 0; // automatic choice for analysis analysis - switch(this->parallel_method) { + switch (this->parallel_method) { case SolverMumpsOptions::_fully_distributed: - icntl(18) = 3; //fully distributed + icntl(18) = 3; // fully distributed - this->mumps_data.nz_loc = matrix->getNbNonZero(); - this->mumps_data.irn_loc = matrix->getIRN().storage(); - this->mumps_data.jcn_loc = matrix->getJCN().storage(); - break; + this->mumps_data.nz_loc = matrix->getNbNonZero(); + this->mumps_data.irn_loc = matrix->getIRN().storage(); + this->mumps_data.jcn_loc = matrix->getJCN().storage(); + break; case SolverMumpsOptions::_not_parallel: case SolverMumpsOptions::_master_slave_distributed: - icntl(18) = 0; //centralized + icntl(18) = 0; // centralized - if(prank == 0) { - this->mumps_data.nz = matrix->getNbNonZero(); + if (prank == 0) { + this->mumps_data.nz = matrix->getNbNonZero(); this->mumps_data.irn = matrix->getIRN().storage(); this->mumps_data.jcn = matrix->getJCN().storage(); } else { - this->mumps_data.nz = 0; + this->mumps_data.nz = 0; this->mumps_data.irn = NULL; this->mumps_data.jcn = NULL; } break; default: AKANTU_DEBUG_ERROR("This case should not happen!!"); } } /* -------------------------------------------------------------------------- */ void SolverMumps::initialize(SolverOptions & options) { AKANTU_DEBUG_IN(); - if(SolverMumpsOptions * opt = dynamic_cast(&options)) { + if (SolverMumpsOptions * opt = dynamic_cast(&options)) { this->parallel_method = opt->parallel_method; } this->mumps_data.par = 1; // The host is part of computations - switch(this->parallel_method) { - case SolverMumpsOptions::_not_parallel: break; + switch (this->parallel_method) { + case SolverMumpsOptions::_not_parallel: + break; case SolverMumpsOptions::_master_slave_distributed: this->mumps_data.par = 0; // The host is not part of the computations case SolverMumpsOptions::_fully_distributed: #ifdef AKANTU_USE_MPI - const StaticCommunicatorMPI & mpi_st_comm = dynamic_cast(communicator.getRealStaticCommunicator()); - this->mumps_data.comm_fortran = MPI_Comm_c2f(mpi_st_comm.getMPITypeWrapper().getMPICommunicator()); + const StaticCommunicatorMPI & mpi_st_comm = + dynamic_cast( + communicator.getRealStaticCommunicator()); + this->mumps_data.comm_fortran = + MPI_Comm_c2f(mpi_st_comm.getMPITypeWrapper().getMPICommunicator()); #endif break; } - this->mumps_data.sym = 2 * (matrix->getSparseMatrixType() == _symmetric); this->prank = communicator.whoAmI(); - this->mumps_data.job = _smj_initialize; //initialize + this->mumps_data.job = _smj_initialize; // initialize dmumps_c(&this->mumps_data); this->is_mumps_data_initialized = true; /* ------------------------------------------------------------------------ */ UInt size = matrix->getSize(); - if(prank == 0) { - std::stringstream sstr_rhs; sstr_rhs << id << ":rhs"; + if (prank == 0) { + std::stringstream sstr_rhs; + sstr_rhs << id << ":rhs"; this->rhs = &(alloc(sstr_rhs.str(), size, 1, 0.)); } else { this->rhs = NULL; } this->mumps_data.nz_alloc = 0; - this->mumps_data.n = size; + this->mumps_data.n = size; /* ------------------------------------------------------------------------ */ // Output setup - if(AKANTU_DEBUG_TEST(dblTrace)) { + if (AKANTU_DEBUG_TEST(dblTrace)) { icntl(1) = 6; icntl(2) = 2; icntl(3) = 2; icntl(4) = 4; } else { /// No outputs icntl(1) = 6; // error output icntl(2) = 0; // dignostics output icntl(3) = 0; // informations icntl(4) = 0; // no outputs } - if(AKANTU_DEBUG_TEST(dblDump)) { + if (AKANTU_DEBUG_TEST(dblDump)) { strcpy(this->mumps_data.write_problem, "mumps_matrix.mtx"); } this->analysis(); -// icntl(14) = 80; + // icntl(14) = 80; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::setRHS(Array & rhs) { - if(prank == 0) { + if (prank == 0) { - //std::copy(rhs.storage(), rhs.storage() + this->rhs->getSize(), this->rhs->storage()); + // std::copy(rhs.storage(), rhs.storage() + this->rhs->getSize(), + // this->rhs->storage()); DebugLevel dbl = debug::getDebugLevel(); debug::setDebugLevel(dblError); matrix->getDOFSynchronizer().gather(rhs, 0, this->rhs); debug::setDebugLevel(dbl); } else { this->matrix->getDOFSynchronizer().gather(rhs, 0); } } /* -------------------------------------------------------------------------- */ void SolverMumps::analysis() { AKANTU_DEBUG_IN(); initMumpsData(); - this->mumps_data.job = _smj_analyze; //analyze + this->mumps_data.job = _smj_analyze; // analyze dmumps_c(&this->mumps_data); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::factorize() { AKANTU_DEBUG_IN(); - if(parallel_method == SolverMumpsOptions::_fully_distributed) - this->mumps_data.a_loc = this->matrix->getA().storage(); + if (parallel_method == SolverMumpsOptions::_fully_distributed) + this->mumps_data.a_loc = this->matrix->getA().storage(); else { - if(prank == 0) - this->mumps_data.a = this->matrix->getA().storage(); + if (prank == 0) + this->mumps_data.a = this->matrix->getA().storage(); } - if(prank == 0) { + if (prank == 0) { this->mumps_data.rhs = this->rhs->storage(); } - this->mumps_data.job = _smj_factorize; // factorize dmumps_c(&this->mumps_data); this->printError(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::solve() { AKANTU_DEBUG_IN(); - if(prank == 0) { + if (prank == 0) { this->mumps_data.rhs = this->rhs->storage(); } this->mumps_data.job = _smj_solve; // solve dmumps_c(&this->mumps_data); this->printError(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::solve(Array & solution) { AKANTU_DEBUG_IN(); this->solve(); - if(prank == 0) { + if (prank == 0) { matrix->getDOFSynchronizer().scatter(solution, 0, this->rhs); -// std::copy(this->rhs->storage(), this->rhs->storage() + this->rhs->getSize(), solution.storage()); + // std::copy(this->rhs->storage(), this->rhs->storage() + + // this->rhs->getSize(), solution.storage()); } else { this->matrix->getDOFSynchronizer().scatter(solution, 0); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::printError() { Int _info_v[2]; - _info_v[0] = info(1); // to get errors + _info_v[0] = info(1); // to get errors _info_v[1] = -info(1); // to get warnings communicator.allReduce(_info_v, 2, _so_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: { + 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 of 10%"); - this->analysis(); - this->factorize(); - this->solve(); + 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 of 10%"); + this->analysis(); + this->factorize(); + this->solve(); } else { - AKANTU_DEBUG_ERROR("The MUMPS workarray is too small INFO(2)=" << info(2) << "No further increase possible"); break; + 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) = " + 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_DEBUG_WARNING("Warning in mumps during solve process, check mumps " + "user guide INFO(1) = " + << _info_v[1]); } } - __END_AKANTU__ diff --git a/src/solver/solver_mumps.hh b/src/solver/solver_mumps.hh index fdbeeadd2..3b422bda0 100644 --- a/src/solver/solver_mumps.hh +++ b/src/solver/solver_mumps.hh @@ -1,155 +1,145 @@ /** * @file solver_mumps.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Oct 16 2015 * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_MUMPS_HH__ #define __AKANTU_SOLVER_MUMPS_HH__ #include #include "solver.hh" #include "static_communicator.hh" __BEGIN_AKANTU__ class SolverMumpsOptions : public SolverOptions { public: enum ParallelMethod { _not_parallel, _fully_distributed, _master_slave_distributed }; - SolverMumpsOptions(ParallelMethod parallel_method = _fully_distributed) : - SolverOptions(), - parallel_method(parallel_method) { } + SolverMumpsOptions(ParallelMethod parallel_method = _fully_distributed) + : SolverOptions(), parallel_method(parallel_method) {} private: friend class SolverMumps; ParallelMethod parallel_method; }; class SolverMumps : public Solver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - - SolverMumps(SparseMatrix & sparse_matrix, - const ID & id = "solver_mumps", - const MemoryID & memory_id = 0); + SolverMumps(SparseMatrix & sparse_matrix, const ID & id = "solver_mumps", + const MemoryID & memory_id = 0); virtual ~SolverMumps(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// build the profile and do the analysis part virtual void initialize(SolverOptions & options = _solver_no_options); void initializeSlave(SolverOptions & options = _solver_no_options); /// analysis (symbolic facto + permutations) virtual void analysis(); /// factorize the matrix virtual void factorize(); /// solve the system virtual void solve(Array & solution); virtual void solve(); void solveSlave(); virtual void setRHS(Array & rhs); private: /// print the error if any happened in mumps void printError(); /// clean the mumps info virtual void destroyInternalData(); /// set internal values; void initMumpsData(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ private: /// access the control variable - inline Int & icntl(UInt i) { - return mumps_data.icntl[i - 1]; - } + 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]; - } + inline Int & info(UInt i) { return mumps_data.info[i - 1]; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: - /// mumps data DMUMPS_STRUC_C mumps_data; /// specify if the mumps_data are initialized or not bool is_mumps_data_initialized; UInt prank; /* ------------------------------------------------------------------------ */ /* Local types */ /* ------------------------------------------------------------------------ */ private: SolverMumpsOptions::ParallelMethod 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 }; }; - __END_AKANTU__ #endif /* __AKANTU_SOLVER_MUMPS_HH__ */ diff --git a/src/solver/solver_petsc.cc b/src/solver/solver_petsc.cc index bd5870127..13996502a 100644 --- a/src/solver/solver_petsc.cc +++ b/src/solver/solver_petsc.cc @@ -1,1110 +1,1170 @@ /** * @file solver_petsc.cc * * @author Alejandro M. Aragón * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Tue May 13 2014 * @date last modification: Fri Oct 16 2015 * * @brief Solver class implementation for the petsc solver * * @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 "solver_petsc.hh" #include "petsc_wrapper.hh" #include "petsc_matrix.hh" #include "static_communicator.hh" #include "static_communicator_mpi.hh" #include "mpi_type_wrapper.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -SolverPETSc::SolverPETSc(SparseMatrix & matrix, - const ID & id, - const MemoryID & memory_id) : - Solver(matrix, id, memory_id), is_petsc_data_initialized(false), - petsc_solver_wrapper(new PETScSolverWrapper) { -} +SolverPETSc::SolverPETSc(SparseMatrix & matrix, const ID & id, + const MemoryID & memory_id) + : Solver(matrix, id, memory_id), is_petsc_data_initialized(false), + petsc_solver_wrapper(new PETScSolverWrapper) {} /* -------------------------------------------------------------------------- */ SolverPETSc::~SolverPETSc() { AKANTU_DEBUG_IN(); this->destroyInternalData(); AKANTU_DEBUG_OUT(); - } +} /* -------------------------------------------------------------------------- */ void SolverPETSc::destroyInternalData() { AKANTU_DEBUG_IN(); - if(this->is_petsc_data_initialized) { + if (this->is_petsc_data_initialized) { PetscErrorCode ierr; - ierr = KSPDestroy(&(this->petsc_solver_wrapper->ksp)); CHKERRXX(ierr); - ierr = VecDestroy(&(this->petsc_solver_wrapper->rhs)); CHKERRXX(ierr); - ierr = VecDestroy(&(this->petsc_solver_wrapper->solution)); CHKERRXX(ierr); + ierr = KSPDestroy(&(this->petsc_solver_wrapper->ksp)); + CHKERRXX(ierr); + ierr = VecDestroy(&(this->petsc_solver_wrapper->rhs)); + CHKERRXX(ierr); + ierr = VecDestroy(&(this->petsc_solver_wrapper->solution)); + CHKERRXX(ierr); delete petsc_solver_wrapper; this->is_petsc_data_initialized = false; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverPETSc::initialize(SolverOptions & options) { AKANTU_DEBUG_IN(); #if defined(AKANTU_USE_MPI) - StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); - const StaticCommunicatorMPI & mpi_st_comm = - dynamic_cast(comm.getRealStaticCommunicator()); - this->petsc_solver_wrapper->communicator = mpi_st_comm.getMPITypeWrapper().getMPICommunicator(); + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + const StaticCommunicatorMPI & mpi_st_comm = + dynamic_cast( + comm.getRealStaticCommunicator()); + this->petsc_solver_wrapper->communicator = + mpi_st_comm.getMPITypeWrapper().getMPICommunicator(); #else - this->petsc_solver_wrapper->communicator = PETSC_COMM_SELF; + this->petsc_solver_wrapper->communicator = PETSC_COMM_SELF; #endif PetscErrorCode ierr; - PETScMatrix * petsc_matrix = static_cast(this->matrix); + PETScMatrix * petsc_matrix = static_cast(this->matrix); /// create a solver context - ierr = KSPCreate(this->petsc_solver_wrapper->communicator, &(this->petsc_solver_wrapper->ksp)); CHKERRXX(ierr); - + ierr = KSPCreate(this->petsc_solver_wrapper->communicator, + &(this->petsc_solver_wrapper->ksp)); + CHKERRXX(ierr); + /// create the PETSc vector for the right-hand side - ierr = VecCreate(this->petsc_solver_wrapper->communicator, &(this->petsc_solver_wrapper->rhs)); CHKERRXX(ierr); - + ierr = VecCreate(this->petsc_solver_wrapper->communicator, + &(this->petsc_solver_wrapper->rhs)); + CHKERRXX(ierr); + ierr = VecSetSizes((this->petsc_solver_wrapper->rhs), - petsc_matrix->getLocalSize(), - petsc_matrix->getSize()); CHKERRXX(ierr); - ierr = VecSetFromOptions((this->petsc_solver_wrapper->rhs)); CHKERRXX(ierr); + petsc_matrix->getLocalSize(), petsc_matrix->getSize()); + CHKERRXX(ierr); + ierr = VecSetFromOptions((this->petsc_solver_wrapper->rhs)); + CHKERRXX(ierr); /// create the PETSc vector for the solution - ierr = VecDuplicate((this->petsc_solver_wrapper->rhs), &(this->petsc_solver_wrapper->solution)); CHKERRXX(ierr); + ierr = VecDuplicate((this->petsc_solver_wrapper->rhs), + &(this->petsc_solver_wrapper->solution)); + CHKERRXX(ierr); /// set the solution to zero ierr = VecZeroEntries(this->petsc_solver_wrapper->solution); - this->is_petsc_data_initialized = true; + this->is_petsc_data_initialized = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverPETSc::setRHS(Array & rhs) { PetscErrorCode ierr; - PETScMatrix * petsc_matrix = static_cast(this->matrix); + PETScMatrix * petsc_matrix = static_cast(this->matrix); UInt nb_component = rhs.getNbComponent(); UInt size = rhs.getSize(); for (UInt i = 0; i < size; ++i) { for (UInt j = 0; j < nb_component; ++j) { UInt i_local = i * nb_component + j; if (this->matrix->getDOFSynchronizer().isLocalOrMasterDOF(i_local)) { - Int i_global = this->matrix->getDOFSynchronizer().getDOFGlobalID(i_local); - AOApplicationToPetsc(petsc_matrix->getPETScMatrixWrapper()->ao, 1, &(i_global) ); - ierr = VecSetValue((this->petsc_solver_wrapper->rhs), i_global, rhs(i,j), INSERT_VALUES); CHKERRXX(ierr); + Int i_global = + this->matrix->getDOFSynchronizer().getDOFGlobalID(i_local); + AOApplicationToPetsc(petsc_matrix->getPETScMatrixWrapper()->ao, 1, + &(i_global)); + ierr = VecSetValue((this->petsc_solver_wrapper->rhs), i_global, + rhs(i, j), INSERT_VALUES); + CHKERRXX(ierr); } } } - ierr = VecAssemblyBegin(this->petsc_solver_wrapper->rhs); CHKERRXX(ierr); - ierr = VecAssemblyEnd(this->petsc_solver_wrapper->rhs); CHKERRXX(ierr); - /// ierr = VecCopy((this->petsc_solver_wrapper->rhs), (this->petsc_solver_wrapper->solution)); CHKERRXX(ierr); - - + ierr = VecAssemblyBegin(this->petsc_solver_wrapper->rhs); + CHKERRXX(ierr); + ierr = VecAssemblyEnd(this->petsc_solver_wrapper->rhs); + CHKERRXX(ierr); + /// ierr = VecCopy((this->petsc_solver_wrapper->rhs), + /// (this->petsc_solver_wrapper->solution)); CHKERRXX(ierr); } /* -------------------------------------------------------------------------- */ void SolverPETSc::solve() { AKANTU_DEBUG_IN(); PetscErrorCode ierr; - ierr = KSPSolve(this->petsc_solver_wrapper->ksp, this->petsc_solver_wrapper->rhs, this->petsc_solver_wrapper->solution); CHKERRXX(ierr); + ierr = + KSPSolve(this->petsc_solver_wrapper->ksp, this->petsc_solver_wrapper->rhs, + this->petsc_solver_wrapper->solution); + CHKERRXX(ierr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverPETSc::solve(Array & solution) { AKANTU_DEBUG_IN(); this->solution = &solution; this->solve(); PetscErrorCode ierr; - PETScMatrix * petsc_matrix = static_cast(this->matrix); + PETScMatrix * petsc_matrix = static_cast(this->matrix); - - // ierr = VecGetOwnershipRange(this->petsc_solver_wrapper->solution, solution_begin, solution_end); CHKERRXX(ierr); - // ierr = VecGetArray(this->petsc_solver_wrapper->solution, PetscScalar **array); CHKERRXX(ierr); + // ierr = VecGetOwnershipRange(this->petsc_solver_wrapper->solution, + // solution_begin, solution_end); CHKERRXX(ierr); + // ierr = VecGetArray(this->petsc_solver_wrapper->solution, PetscScalar + // **array); CHKERRXX(ierr); UInt nb_component = solution.getNbComponent(); UInt size = solution.getSize(); - + for (UInt i = 0; i < size; ++i) { for (UInt j = 0; j < nb_component; ++j) { UInt i_local = i * nb_component + j; if (this->matrix->getDOFSynchronizer().isLocalOrMasterDOF(i_local)) { - Int i_global = this->matrix->getDOFSynchronizer().getDOFGlobalID(i_local); - ierr = AOApplicationToPetsc(petsc_matrix->getPETScMatrixWrapper()->ao, 1, &(i_global) ); CHKERRXX(ierr); - ierr = VecGetValues(this->petsc_solver_wrapper->solution, 1, &i_global, &solution(i,j)); CHKERRXX(ierr); - + Int i_global = + this->matrix->getDOFSynchronizer().getDOFGlobalID(i_local); + ierr = AOApplicationToPetsc(petsc_matrix->getPETScMatrixWrapper()->ao, + 1, &(i_global)); + CHKERRXX(ierr); + ierr = VecGetValues(this->petsc_solver_wrapper->solution, 1, &i_global, + &solution(i, j)); + CHKERRXX(ierr); } } } synch_registry->synchronize(_gst_solver_solution); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverPETSc::setOperators() { AKANTU_DEBUG_IN(); PetscErrorCode ierr; - /// set the matrix that defines the linear system and the matrix for preconditioning (here they are the same) - PETScMatrix * petsc_matrix = static_cast(this->matrix); + /// set the matrix that defines the linear system and the matrix for + /// preconditioning (here they are the same) + PETScMatrix * petsc_matrix = static_cast(this->matrix); -#if PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 5 - ierr = KSPSetOperators(this->petsc_solver_wrapper->ksp, petsc_matrix->getPETScMatrixWrapper()->mat, petsc_matrix->getPETScMatrixWrapper()->mat); CHKERRXX(ierr); +#if PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 5 + ierr = KSPSetOperators(this->petsc_solver_wrapper->ksp, + petsc_matrix->getPETScMatrixWrapper()->mat, + petsc_matrix->getPETScMatrixWrapper()->mat); + CHKERRXX(ierr); #else - ierr = KSPSetOperators(this->petsc_solver_wrapper->ksp, petsc_matrix->getPETScMatrixWrapper()->mat, petsc_matrix->getPETScMatrixWrapper()->mat, SAME_NONZERO_PATTERN); CHKERRXX(ierr); + ierr = KSPSetOperators(this->petsc_solver_wrapper->ksp, + petsc_matrix->getPETScMatrixWrapper()->mat, + petsc_matrix->getPETScMatrixWrapper()->mat, + SAME_NONZERO_PATTERN); + CHKERRXX(ierr); #endif - /// If this is not called the solution vector is zeroed in the call to KSPSolve(). - ierr = KSPSetInitialGuessNonzero(this->petsc_solver_wrapper->ksp, PETSC_TRUE); CHKERRXX(ierr); - ierr = KSPSetFromOptions(this->petsc_solver_wrapper->ksp); CHKERRXX(ierr); + /// If this is not called the solution vector is zeroed in the call to + /// KSPSolve(). + ierr = KSPSetInitialGuessNonzero(this->petsc_solver_wrapper->ksp, PETSC_TRUE); + CHKERRXX(ierr); + ierr = KSPSetFromOptions(this->petsc_solver_wrapper->ksp); + CHKERRXX(ierr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void finalize_petsc() { - + // static bool finalized = false; // if (!finalized) { - + // cout<<"*** INFO *** PETSc is finalizing..."<first; // ierr = VecSetValues(b, 1, &row, &it->second, ADD_VALUES); // } - + // #ifdef CPPUTILS_VERBOSE // out<<" Right hand side set... "< subs = AA.unhash(it->first); // PetscInt row = subs.first; // PetscInt col = subs.second; -// ierr = MatSetValues(A_, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERRCONTINUE(ierr); +// ierr = MatSetValues(A_, 1, &row, 1, &col, &it->second, +// ADD_VALUES);CHKERRCONTINUE(ierr); // } - + // #ifdef CPPUTILS_VERBOSE // out<<" Filled sparse matrix... "< -pc_type -ksp_monitor -ksp_rtol // // These options will override those specified above as long as // // KSPSetFromOptions() is called _after_ any other customization // // routines. // // */ // // ierr = KSPSetFromOptions(ksp);CHKERRCONTINUE(ierr); - - + // #ifdef CPPUTILS_VERBOSE // out<<" Solving system... "<first; // ierr = VecGetValues(x_all, 1, &row, &val); // if (val != zero) // xx[row] = val; // } - - + // #ifdef CPPUTILS_VERBOSE // out<<" Solution vector copied... "< Kpf_nz = KKpf.non_zero_off_diagonal(); // std::pair Kpp_nz = KKpp.non_zero_off_diagonal(); - -// ierr = MatMPIAIJSetPreallocation(Kpf, Kpf_nz.first, PETSC_NULL, Kpf_nz.second, PETSC_NULL); CHKERRCONTINUE(ierr); -// ierr = MatMPIAIJSetPreallocation(Kpp, Kpp_nz.first, PETSC_NULL, Kpp_nz.second, PETSC_NULL); CHKERRCONTINUE(ierr); -// ierr = MatSeqAIJSetPreallocation(Kpf, Kpf_nz.first, PETSC_NULL); CHKERRCONTINUE(ierr); -// ierr = MatSeqAIJSetPreallocation(Kpp, Kpp_nz.first, PETSC_NULL); CHKERRCONTINUE(ierr); - -// for (sparse_matrix_type::const_hash_iterator it = KKpf.map_.begin(); it != KKpf.map_.end(); ++it) { - + +// ierr = MatMPIAIJSetPreallocation(Kpf, Kpf_nz.first, PETSC_NULL, +// Kpf_nz.second, PETSC_NULL); CHKERRCONTINUE(ierr); +// ierr = MatMPIAIJSetPreallocation(Kpp, Kpp_nz.first, PETSC_NULL, +// Kpp_nz.second, PETSC_NULL); CHKERRCONTINUE(ierr); +// ierr = MatSeqAIJSetPreallocation(Kpf, Kpf_nz.first, PETSC_NULL); +// CHKERRCONTINUE(ierr); +// ierr = MatSeqAIJSetPreallocation(Kpp, Kpp_nz.first, PETSC_NULL); +// CHKERRCONTINUE(ierr); + +// for (sparse_matrix_type::const_hash_iterator it = KKpf.map_.begin(); it != +// KKpf.map_.end(); ++it) { + // // get subscripts // std::pair subs = KKpf.unhash(it->first); // int row = subs.first; // int col = subs.second; -// ierr = MatSetValues(Kpf, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERRCONTINUE(ierr); +// ierr = MatSetValues(Kpf, 1, &row, 1, &col, &it->second, +// ADD_VALUES);CHKERRCONTINUE(ierr); // } - -// for (sparse_matrix_type::const_hash_iterator it = KKpp.map_.begin(); it != KKpp.map_.end(); ++it) { - + +// for (sparse_matrix_type::const_hash_iterator it = KKpp.map_.begin(); it != +// KKpp.map_.end(); ++it) { + // // get subscripts // std::pair subs = KKpp.unhash(it->first); // int row = subs.first; // int col = subs.second; -// ierr = MatSetValues(Kpf, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERRCONTINUE(ierr); +// ierr = MatSetValues(Kpf, 1, &row, 1, &col, &it->second, +// ADD_VALUES);CHKERRCONTINUE(ierr); // } - + // #ifdef CPPUTILS_VERBOSE // out<<" Filled sparse matrices... "<first; // ierr = VecSetValues(Up, 1, &row, &it->second, ADD_VALUES); // } - + // /* // Assemble vector, using the 2-step process: // VecAssemblyBegin(), VecAssemblyEnd() // Computations can be done while messages are in transition // by placing code between these two statements. // */ // ierr = VecAssemblyBegin(Up);CHKERRCONTINUE(ierr); // ierr = VecAssemblyEnd(Up);CHKERRCONTINUE(ierr); - + // // add Kpf*Uf // ierr = MatMult(Kpf, x_, Pf); - + // // add Kpp*Up // ierr = MatMultAdd(Kpp, Up, Pf, Pp); - + // #ifdef CPPUTILS_VERBOSE // out<<" Matrices multiplied... "<first; // ierr = VecSetValues(r, 1, &row, &it->second, ADD_VALUES); // } -// for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) { +// for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != +// bb.map_.end(); ++it) { // int row = it->first; // ierr = VecSetValues(r, 1, &row, &it->second, ADD_VALUES); // } - + // /* // Assemble vector, using the 2-step process: // VecAssemblyBegin(), VecAssemblyEnd() // Computations can be done while messages are in transition // by placing code between these two statements. // */ // ierr = VecAssemblyBegin(r);CHKERRCONTINUE(ierr); // ierr = VecAssemblyEnd(r);CHKERRCONTINUE(ierr); - - + // sparse_vector_type rr(aa.size()); - -// for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != aa.map_.end(); ++it) { + +// for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != +// aa.map_.end(); ++it) { // int row = it->first; // ierr = VecGetValues(r, 1, &row, &rr[row]); // } -// for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) { +// for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != +// bb.map_.end(); ++it) { // int row = it->first; // ierr = VecGetValues(r, 1, &row, &rr[row]); // } - + // #ifdef CPPUTILS_VERBOSE // out<<" Vector copied... "< subs = aa.unhash(it->first); // int row = subs.first; // int col = subs.second; - + // if (flag == Add_t) // ierr = MatSetValues(r, 1, &row, 1, &col, &it->second, ADD_VALUES); // else if (flag == Insert_t) // ierr = MatSetValues(r, 1, &row, 1, &col, &it->second, INSERT_VALUES); // CHKERRCONTINUE(ierr); // } - + // #ifdef CPPUTILS_VERBOSE // out<<" Matrix filled..."<first; // if (flag == Add_t) // ierr = VecSetValues(r, 1, &row, &it->second, ADD_VALUES); // else if (flag == Insert_t) // ierr = VecSetValues(r, 1, &row, &it->second, INSERT_VALUES); // CHKERRCONTINUE(ierr); // } - + // #ifdef CPPUTILS_VERBOSE // out<<" Vector filled..."<(comm.getRealStaticCommunicator()); -// // if(mumps_data.comm_fortran == MPI_Comm_c2f(comm_mpi.getMPICommunicator())) +// // dynamic_cast(comm.getRealStaticCommunicator()); +// // if(mumps_data.comm_fortran == +// MPI_Comm_c2f(comm_mpi.getMPICommunicator())) // // destroyMumpsData(); // // } catch(...) {} -// // +// // // // AKANTU_DEBUG_OUT(); // //} // // -// ///* -------------------------------------------------------------------------- */ -// //void SolverMumps::initMumpsData(SolverMumpsOptions::ParallelMethod parallel_method) { +// ///* +// -------------------------------------------------------------------------- */ +// //void SolverMumps::initMumpsData(SolverMumpsOptions::ParallelMethod +// parallel_method) { // // switch(parallel_method) { // // case SolverMumpsOptions::_fully_distributed: // // icntl(18) = 3; //fully distributed // // icntl(28) = 0; //automatic choice -// // +// // // // mumps_data.nz_loc = matrix->getNbNonZero(); // // mumps_data.irn_loc = matrix->getIRN().values; // // mumps_data.jcn_loc = matrix->getJCN().values; // // break; // // case SolverMumpsOptions::_master_slave_distributed: // // if(prank == 0) { // // mumps_data.nz = matrix->getNbNonZero(); // // mumps_data.irn = matrix->getIRN().values; // // mumps_data.jcn = matrix->getJCN().values; // // } else { // // mumps_data.nz = 0; // // mumps_data.irn = NULL; // // mumps_data.jcn = NULL; -// // +// // // // icntl(18) = 0; //centralized // // icntl(28) = 0; //sequential analysis // // } // // break; // // } // //} // // -// ///* -------------------------------------------------------------------------- */ +// ///* +// -------------------------------------------------------------------------- */ // //void SolverMumps::initialize(SolverOptions & options) { // // AKANTU_DEBUG_IN(); -// // +// // // // mumps_data.par = 1; -// // -// // if(SolverMumpsOptions * opt = dynamic_cast(&options)) { -// // if(opt->parallel_method == SolverMumpsOptions::_master_slave_distributed) { +// // +// // if(SolverMumpsOptions * opt = dynamic_cast(&options)) { +// // if(opt->parallel_method == +// SolverMumpsOptions::_master_slave_distributed) { // // mumps_data.par = 0; // // } // // } -// // +// // // // mumps_data.sym = 2 * (matrix->getSparseMatrixType() == _symmetric); // // prank = communicator.whoAmI(); // //#ifdef AKANTU_USE_MPI -// // mumps_data.comm_fortran = MPI_Comm_c2f(dynamic_cast(communicator.getRealStaticCommunicator()).getMPICommunicator()); +// // mumps_data.comm_fortran = MPI_Comm_c2f(dynamic_cast(communicator.getRealStaticCommunicator()).getMPICommunicator()); // //#endif -// // +// // // // if(AKANTU_DEBUG_TEST(dblTrace)) { // // icntl(1) = 2; // // icntl(2) = 2; // // icntl(3) = 2; // // icntl(4) = 4; // // } -// // +// // // // mumps_data.job = _smj_initialize; //initialize // // dmumps_c(&mumps_data); // // is_mumps_data_initialized = true; -// // -// // /* ------------------------------------------------------------------------ */ +// // +// // /* +// ------------------------------------------------------------------------ */ // // UInt size = matrix->getSize(); -// // +// // // // if(prank == 0) { // // std::stringstream sstr_rhs; sstr_rhs << id << ":rhs"; // // rhs = &(alloc(sstr_rhs.str(), size, 1, REAL_INIT_VALUE)); // // } else { // // rhs = NULL; // // } -// // +// // // // /// No outputs // // icntl(1) = 0; // // icntl(2) = 0; // // icntl(3) = 0; // // icntl(4) = 0; // // mumps_data.nz_alloc = 0; -// // +// // // // if (AKANTU_DEBUG_TEST(dblDump)) icntl(4) = 4; -// // +// // // // mumps_data.n = size; -// // +// // // // if(AKANTU_DEBUG_TEST(dblDump)) { // // strcpy(mumps_data.write_problem, "mumps_matrix.mtx"); // // } -// // -// // /* ------------------------------------------------------------------------ */ +// // +// // /* +// ------------------------------------------------------------------------ */ // // // Default Scaling // // icntl(8) = 77; -// // +// // // // icntl(5) = 0; // Assembled matrix -// // +// // // // SolverMumpsOptions * opt = dynamic_cast(&options); // // if(opt) // // parallel_method = opt->parallel_method; -// // +// // // // initMumpsData(parallel_method); -// // +// // // // mumps_data.job = _smj_analyze; //analyze // // dmumps_c(&mumps_data); -// // +// // // // AKANTU_DEBUG_OUT(); // //} // // -// ///* -------------------------------------------------------------------------- */ +// ///* +// -------------------------------------------------------------------------- */ // //void SolverMumps::setRHS(Array & rhs) { // // if(prank == 0) { // // matrix->getDOFSynchronizer().gather(rhs, 0, this->rhs); // // } else { // // matrix->getDOFSynchronizer().gather(rhs, 0); // // } // //} // // -// ///* -------------------------------------------------------------------------- */ +// ///* +// -------------------------------------------------------------------------- */ // //void SolverMumps::solve() { // // AKANTU_DEBUG_IN(); -// // +// // // // if(parallel_method == SolverMumpsOptions::_fully_distributed) // // mumps_data.a_loc = matrix->getA().values; // // else // // if(prank == 0) { // // mumps_data.a = matrix->getA().values; // // } -// // +// // // // if(prank == 0) { // // mumps_data.rhs = rhs->values; // // } -// // +// // // // /// Default centralized dense second member // // icntl(20) = 0; // // icntl(21) = 0; -// // +// // // // mumps_data.job = _smj_factorize_solve; //solve // // dmumps_c(&mumps_data); -// // +// // // // AKANTU_DEBUG_ASSERT(info(1) != -10, "Singular matrix"); // // AKANTU_DEBUG_ASSERT(info(1) == 0, -// // "Error in mumps during solve process, check mumps user guide INFO(1) =" +// // "Error in mumps during solve process, check mumps +// user guide INFO(1) =" // // << info(1)); -// // +// // // // AKANTU_DEBUG_OUT(); // //} // // -// ///* -------------------------------------------------------------------------- */ +// ///* +// -------------------------------------------------------------------------- */ // //void SolverMumps::solve(Array & solution) { // // AKANTU_DEBUG_IN(); -// // +// // // // solve(); -// // +// // // // if(prank == 0) { // // matrix->getDOFSynchronizer().scatter(solution, 0, this->rhs); // // } else { // // matrix->getDOFSynchronizer().scatter(solution, 0); // // } -// // +// // // // AKANTU_DEBUG_OUT(); // //} __END_AKANTU__ - diff --git a/src/solver/static_solver.cc b/src/solver/static_solver.cc index 6f5c43a9c..150bf6080 100644 --- a/src/solver/static_solver.cc +++ b/src/solver/static_solver.cc @@ -1,116 +1,123 @@ /** * @file static_solver.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Mon Oct 13 2014 * @date last modification: Fri Oct 16 2015 * * @brief implementation of the static solver * * @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 "static_solver.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_PETSC #include #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -StaticSolver::StaticSolver() : CommunicatorEventHandler(), is_initialized(false) { +StaticSolver::StaticSolver() + : CommunicatorEventHandler(), is_initialized(false) { StaticCommunicator::getStaticCommunicator().registerEventHandler(*this); } - /* -------------------------------------------------------------------------- */ StaticSolver::~StaticSolver() { --this->nb_references; - if(this->nb_references == 0) { + if (this->nb_references == 0) { StaticCommunicator::getStaticCommunicator().unregisterEventHandler(*this); delete this->static_solver; } } /* -------------------------------------------------------------------------- */ StaticSolver & StaticSolver::getStaticSolver() { - if(nb_references == 0) + if (nb_references == 0) static_solver = new StaticSolver(); ++nb_references; return *static_solver; } #ifdef AKANTU_USE_PETSC #if PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 5 -static PetscErrorCode PETScErrorHandler(MPI_Comm, - int line, const char * dir, const char *file, +static PetscErrorCode PETScErrorHandler(MPI_Comm, int line, const char * dir, + const char * file, PetscErrorCode number, PetscErrorType type, - const char *message, - void *) { - AKANTU_DEBUG_ERROR("An error occured in PETSc in file \"" << file << ":" << line << "\" - PetscErrorCode "<< number << " - \""<< message << "\""); + const char * message, void *) { + AKANTU_DEBUG_ERROR("An error occured in PETSc in file \"" + << file << ":" << line << "\" - PetscErrorCode " << number + << " - \"" << message << "\""); } #else -static PetscErrorCode PETScErrorHandler(MPI_Comm, - int line, const char * func, const char * dir, const char *file, +static PetscErrorCode PETScErrorHandler(MPI_Comm, int line, const char * func, + const char * dir, const char * file, PetscErrorCode number, PetscErrorType type, - const char *message, - void *) { - AKANTU_DEBUG_ERROR("An error occured in PETSc in file \"" << file << ":" << line << "\" - PetscErrorCode "<< number << " - \""<< message << "\""); + const char * message, void *) { + AKANTU_DEBUG_ERROR("An error occured in PETSc in file \"" + << file << ":" << line << "\" - PetscErrorCode " << number + << " - \"" << message << "\""); } #endif #endif /* -------------------------------------------------------------------------- */ -void StaticSolver::initialize(int & argc, char ** & argv) { - if (this->is_initialized) return; - // AKANTU_DEBUG_ASSERT(this->is_initialized != true, "The static solver has already been initialized"); +void StaticSolver::initialize(int & argc, char **& argv) { + if (this->is_initialized) + return; +// AKANTU_DEBUG_ASSERT(this->is_initialized != true, "The static solver has +// already been initialized"); #ifdef AKANTU_USE_PETSC - PetscErrorCode petsc_error = PetscInitialize(&argc, &argv, NULL, NULL); - if(petsc_error != 0) { - AKANTU_DEBUG_ERROR("An error occured while initializing Petsc (PetscErrorCode "<< petsc_error << ")"); - } - PetscPushErrorHandler(PETScErrorHandler, NULL); + PetscErrorCode petsc_error = PetscInitialize(&argc, &argv, NULL, NULL); + if (petsc_error != 0) { + AKANTU_DEBUG_ERROR( + "An error occured while initializing Petsc (PetscErrorCode " + << petsc_error << ")"); + } + PetscPushErrorHandler(PETScErrorHandler, NULL); #endif - this->is_initialized = true; - } + this->is_initialized = true; +} /* -------------------------------------------------------------------------- */ void StaticSolver::finalize() { - ParentEventHandler::sendEvent(StaticSolverEvent::BeforeStaticSolverDestroyEvent()); - + ParentEventHandler::sendEvent( + StaticSolverEvent::BeforeStaticSolverDestroyEvent()); - AKANTU_DEBUG_ASSERT(this->is_initialized == true, "The static solver has not been initialized"); + AKANTU_DEBUG_ASSERT(this->is_initialized == true, + "The static solver has not been initialized"); #ifdef AKANTU_USE_PETSC PetscFinalize(); #endif this->is_initialized = false; } __END_AKANTU__ diff --git a/src/synchronizer/static_communicator.cc b/src/synchronizer/static_communicator.cc index cd0858d0b..2d040274f 100644 --- a/src/synchronizer/static_communicator.cc +++ b/src/synchronizer/static_communicator.cc @@ -1,116 +1,125 @@ /** * @file static_communicator.cc * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Jan 13 2016 * * @brief implementation of the common part of the static communicator * * @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 "static_communicator.hh" #include "static_communicator_dummy.hh" #ifdef AKANTU_USE_MPI -# include "static_communicator_mpi.hh" +#include "static_communicator_mpi.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ bool StaticCommunicator::is_instantiated = false; StaticCommunicator * StaticCommunicator::static_communicator = NULL; UInt CommunicationRequest::counter = 0; /* -------------------------------------------------------------------------- */ -CommunicationRequest::CommunicationRequest(UInt source, UInt dest) : - source(source), destination(dest) { +CommunicationRequest::CommunicationRequest(UInt source, UInt dest) + : source(source), destination(dest) { this->id = counter++; } /* -------------------------------------------------------------------------- */ -CommunicationRequest::~CommunicationRequest() { - -} +CommunicationRequest::~CommunicationRequest() {} /* -------------------------------------------------------------------------- */ void CommunicationRequest::printself(std::ostream & stream, int indent) const { std::string space; - for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); + for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) + ; stream << space << "CommunicationRequest [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + source : " << source << std::endl; stream << space << " + destination : " << destination << std::endl; stream << space << "]" << std::endl; } - - /* -------------------------------------------------------------------------- */ -StaticCommunicator::StaticCommunicator(int & argc, - char ** & argv, - CommunicatorType type) { +StaticCommunicator::StaticCommunicator(int & argc, char **& argv, + CommunicatorType type) { real_type = type; #ifdef AKANTU_USE_MPI - if(type == _communicator_mpi) { - real_static_communicator = new StaticCommunicatorMPI(argc, argv); - } else { + if (type == _communicator_mpi) { + real_static_communicator = new StaticCommunicatorMPI(argc, argv); + } else { #endif - real_static_communicator = new StaticCommunicatorDummy(argc, argv); + real_static_communicator = new StaticCommunicatorDummy(argc, argv); #ifdef AKANTU_USE_MPI - } + } #endif } /* -------------------------------------------------------------------------- */ -StaticCommunicator & StaticCommunicator::getStaticCommunicator(CommunicatorType type) { +StaticCommunicator::~StaticCommunicator() { + FinalizeCommunicatorEvent * event = new FinalizeCommunicatorEvent(*this); + this->sendEvent(*event); + + this->barrier(); + delete event; + delete real_static_communicator; + is_instantiated = false; + StaticCommunicator::static_communicator = NULL; +} + +/* -------------------------------------------------------------------------- */ +StaticCommunicator & +StaticCommunicator::getStaticCommunicator(CommunicatorType type) { AKANTU_DEBUG_IN(); if (!static_communicator) { int nb_args = 0; char ** null; static_communicator = new StaticCommunicator(nb_args, null, type); + is_instantiated = true; } - is_instantiated = true; AKANTU_DEBUG_OUT(); return *static_communicator; } /* -------------------------------------------------------------------------- */ -StaticCommunicator & StaticCommunicator::getStaticCommunicator(int & argc, - char ** & argv, - CommunicatorType type) { +StaticCommunicator & +StaticCommunicator::getStaticCommunicator(int & argc, char **& argv, + CommunicatorType type) { if (!static_communicator) static_communicator = new StaticCommunicator(argc, argv, type); return getStaticCommunicator(type); } __END_AKANTU__ diff --git a/src/synchronizer/static_communicator.hh b/src/synchronizer/static_communicator.hh index 49e717dfa..9810da691 100644 --- a/src/synchronizer/static_communicator.hh +++ b/src/synchronizer/static_communicator.hh @@ -1,216 +1,216 @@ /** * @file static_communicator.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Jan 13 2016 * * @brief Class handling the parallel communications * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STATIC_COMMUNICATOR_HH__ #define __AKANTU_STATIC_COMMUNICATOR_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_event_handler_manager.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_COMMUNICATOR_LIST_0 BOOST_PP_SEQ_NIL #include "static_communicator_dummy.hh" -#define AKANTU_COMMUNICATOR_LIST_1 \ - BOOST_PP_SEQ_PUSH_BACK(AKANTU_COMMUNICATOR_LIST_0, \ - (_communicator_dummy, (StaticCommunicatorDummy, BOOST_PP_NIL))) +#define AKANTU_COMMUNICATOR_LIST_1 \ + BOOST_PP_SEQ_PUSH_BACK( \ + AKANTU_COMMUNICATOR_LIST_0, \ + (_communicator_dummy, (StaticCommunicatorDummy, BOOST_PP_NIL))) #if defined(AKANTU_USE_MPI) -# include "static_communicator_mpi.hh" -# define AKANTU_COMMUNICATOR_LIST_ALL \ - BOOST_PP_SEQ_PUSH_BACK(AKANTU_COMMUNICATOR_LIST_1, \ - (_communicator_mpi, (StaticCommunicatorMPI, BOOST_PP_NIL))) +#include "static_communicator_mpi.hh" +#define AKANTU_COMMUNICATOR_LIST_ALL \ + BOOST_PP_SEQ_PUSH_BACK( \ + AKANTU_COMMUNICATOR_LIST_1, \ + (_communicator_mpi, (StaticCommunicatorMPI, BOOST_PP_NIL))) #else -# define AKANTU_COMMUNICATOR_LIST_ALL AKANTU_COMMUNICATOR_LIST_1 +#define AKANTU_COMMUNICATOR_LIST_ALL AKANTU_COMMUNICATOR_LIST_1 #endif // AKANTU_COMMUNICATOR_LIST #include "real_static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class RealStaticCommunicator; struct FinalizeCommunicatorEvent { - FinalizeCommunicatorEvent(const StaticCommunicator & comm) : communicator(comm) {} + FinalizeCommunicatorEvent(const StaticCommunicator & comm) + : communicator(comm) {} const StaticCommunicator & communicator; }; class CommunicatorEventHandler { public: virtual ~CommunicatorEventHandler() {} - virtual void onCommunicatorFinalize(__attribute__((unused)) const StaticCommunicator & communicator) { } + virtual void onCommunicatorFinalize(__attribute__((unused)) + const StaticCommunicator & communicator) { + } + private: inline void sendEvent(const FinalizeCommunicatorEvent & event) { onCommunicatorFinalize(event.communicator); } - template - friend class EventHandlerManager; + template friend class EventHandlerManager; }; -class StaticCommunicator : public EventHandlerManager{ +class StaticCommunicator + : public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ protected: - StaticCommunicator(int & argc, char ** & argv, + StaticCommunicator(int & argc, char **& argv, CommunicatorType type = _communicator_mpi); public: - virtual ~StaticCommunicator() { - FinalizeCommunicatorEvent *event = new FinalizeCommunicatorEvent(*this); - this->sendEvent(*event); - - delete event; - delete real_static_communicator; - is_instantiated = false; - StaticCommunicator::static_communicator = NULL; - - }; + virtual ~StaticCommunicator(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /* ------------------------------------------------------------------------ */ /* Point to Point */ /* ------------------------------------------------------------------------ */ - template inline void send(T * buffer, Int size, - Int receiver, Int tag); - template inline void receive(T * buffer, Int size, - Int sender, Int tag); - - template inline CommunicationRequest * asyncSend(T * buffer, - Int size, - Int receiver, - Int tag); - template inline CommunicationRequest * asyncReceive(T * buffer, - Int size, - Int sender, - Int tag); - - template inline void probe(Int sender, Int tag, - CommunicationStatus & status); + template + inline void send(T * buffer, Int size, Int receiver, Int tag); + template + inline void receive(T * buffer, Int size, Int sender, Int tag); + + template + inline CommunicationRequest * asyncSend(T * buffer, Int size, Int receiver, + Int tag); + template + inline CommunicationRequest * asyncReceive(T * buffer, Int size, Int sender, + Int tag); + + template + inline void probe(Int sender, Int tag, CommunicationStatus & status); /* ------------------------------------------------------------------------ */ /* Collectives */ /* ------------------------------------------------------------------------ */ - template inline void reduce(T * values, int nb_values, - const SynchronizerOperation & op, - int root = 0); - template inline void allReduce(T * values, int nb_values, - const SynchronizerOperation & op); - - template inline void allGather(T * values, int nb_values); - template inline void allGatherV(T * values, int * nb_values); - - template inline void gather(T * values, int nb_values, - int root = 0); - template inline void gatherV(T * values, int * nb_values, - int root = 0); - template inline void broadcast(T * values, int nb_values, - int root = 0); + template + inline void reduce(T * values, int nb_values, + const SynchronizerOperation & op, int root = 0); + template + inline void allReduce(T * values, int nb_values, + const SynchronizerOperation & op); + + template inline void allGather(T * values, int nb_values); + template inline void allGatherV(T * values, int * nb_values); + + template + inline void gather(T * values, int nb_values, int root = 0); + template + inline void gatherV(T * values, int * nb_values, int root = 0); + template + inline void broadcast(T * values, int nb_values, int root = 0); inline void barrier(); /* ------------------------------------------------------------------------ */ /* Request handling */ /* ------------------------------------------------------------------------ */ inline bool testRequest(CommunicationRequest * request); inline void wait(CommunicationRequest * request); inline void waitAll(std::vector & requests); inline void freeCommunicationRequest(CommunicationRequest * request); - inline void freeCommunicationRequest(std::vector & requests); + inline void + freeCommunicationRequest(std::vector & requests); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: virtual Int getNbProc() const { return real_static_communicator->psize; }; virtual Int whoAmI() const { return real_static_communicator->prank; }; - AKANTU_GET_MACRO(RealStaticCommunicator, *real_static_communicator, const RealStaticCommunicator &); - AKANTU_GET_MACRO_NOT_CONST(RealStaticCommunicator, *real_static_communicator, RealStaticCommunicator &); + AKANTU_GET_MACRO(RealStaticCommunicator, *real_static_communicator, + const RealStaticCommunicator &); + AKANTU_GET_MACRO_NOT_CONST(RealStaticCommunicator, *real_static_communicator, + RealStaticCommunicator &); - template - Comm & getRealStaticCommunicator() { return dynamic_cast(*real_static_communicator); } + template Comm & getRealStaticCommunicator() { + return dynamic_cast(*real_static_communicator); + } - static StaticCommunicator & getStaticCommunicator(CommunicatorType type = _communicator_mpi); + static StaticCommunicator & + getStaticCommunicator(CommunicatorType type = _communicator_mpi); - static StaticCommunicator & getStaticCommunicator(int & argc, char ** & argv, - CommunicatorType type = _communicator_mpi); + static StaticCommunicator & + getStaticCommunicator(int & argc, char **& argv, + CommunicatorType type = _communicator_mpi); static bool isInstantiated() { return is_instantiated; }; int getMaxTag(); int getMinTag(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: static bool is_instantiated; static StaticCommunicator * static_communicator; RealStaticCommunicator * real_static_communicator; CommunicatorType real_type; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "static_communicator_inline_impl.hh" /* -------------------------------------------------------------------------- */ /* Inline Functions ArrayBase */ /* -------------------------------------------------------------------------- */ -inline std::ostream & operator<<(std::ostream & stream, const CommunicationRequest & _this) -{ +inline std::ostream & operator<<(std::ostream & stream, + const CommunicationRequest & _this) { _this.printself(stream); return stream; } - __END_AKANTU__ #endif /* __AKANTU_STATIC_COMMUNICATOR_HH__ */ diff --git a/src/synchronizer/static_communicator_mpi.cc b/src/synchronizer/static_communicator_mpi.cc index 70c586b3c..a59e8ffac 100644 --- a/src/synchronizer/static_communicator_mpi.cc +++ b/src/synchronizer/static_communicator_mpi.cc @@ -1,496 +1,502 @@ /** * @file static_communicator_mpi.cc * * @author Nicolas Richart * * @date creation: Sun Sep 26 2010 * @date last modification: Wed Jan 13 2016 * * @brief StaticCommunicatorMPI implementation * * @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 "static_communicator_mpi.hh" #include "mpi_type_wrapper.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ MPI_Op MPITypeWrapper::synchronizer_operation_to_mpi_op[_so_null + 1] = { MPI_SUM, MPI_MIN, MPI_MAX, MPI_PROD, MPI_LAND, MPI_BAND, MPI_LOR, MPI_BOR, MPI_LXOR, MPI_BXOR, MPI_MINLOC, MPI_MAXLOC, MPI_OP_NULL}; class CommunicationRequestMPI : public CommunicationRequest { public: CommunicationRequestMPI(UInt source, UInt dest); ~CommunicationRequestMPI(); MPI_Request * getMPIRequest() { return request; }; private: MPI_Request * request; }; /* -------------------------------------------------------------------------- */ /* Implementation */ /* -------------------------------------------------------------------------- */ CommunicationRequestMPI::CommunicationRequestMPI(UInt source, UInt dest) : CommunicationRequest(source, dest) { request = new MPI_Request; } /* -------------------------------------------------------------------------- */ CommunicationRequestMPI::~CommunicationRequestMPI() { delete request; } /* -------------------------------------------------------------------------- */ StaticCommunicatorMPI::StaticCommunicatorMPI(int & argc, char **& argv) : RealStaticCommunicator(argc, argv) { int is_initialized = false; MPI_Initialized(&is_initialized); if (!is_initialized) { MPI_Init(&argc, &argv); } + this->is_externaly_initialized = is_initialized; + mpi_data = new MPITypeWrapper(*this); mpi_data->setMPICommunicator(MPI_COMM_WORLD); } /* -------------------------------------------------------------------------- */ -StaticCommunicatorMPI::~StaticCommunicatorMPI() { MPI_Finalize(); } +StaticCommunicatorMPI::~StaticCommunicatorMPI() { + if (!this->is_externaly_initialized) { + MPI_Finalize(); + } +} /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::send(T * buffer, Int size, Int receiver, Int tag) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Send(buffer, size, type, receiver, tag, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Send."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::receive(T * buffer, Int size, Int sender, Int tag) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Status status; MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Recv(buffer, size, type, sender, tag, communicator, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Recv."); } /* -------------------------------------------------------------------------- */ template CommunicationRequest * StaticCommunicatorMPI::asyncSend(T * buffer, Int size, Int receiver, Int tag) { MPI_Comm communicator = mpi_data->getMPICommunicator(); CommunicationRequestMPI * request = new CommunicationRequestMPI(prank, receiver); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Isend(buffer, size, type, receiver, tag, communicator, request->getMPIRequest()); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Isend."); return request; } /* -------------------------------------------------------------------------- */ template CommunicationRequest * StaticCommunicatorMPI::asyncReceive(T * buffer, Int size, Int sender, Int tag) { MPI_Comm communicator = mpi_data->getMPICommunicator(); CommunicationRequestMPI * request = new CommunicationRequestMPI(sender, prank); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Irecv(buffer, size, type, sender, tag, communicator, request->getMPIRequest()); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Irecv."); return request; } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::probe(Int sender, Int tag, CommunicationStatus & status) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Status mpi_status; #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Probe(sender, tag, communicator, &mpi_status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Probe."); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); int count; MPI_Get_count(&mpi_status, type, &count); status.setSource(mpi_status.MPI_SOURCE); status.setTag(mpi_status.MPI_TAG); status.setSize(count); } /* -------------------------------------------------------------------------- */ bool StaticCommunicatorMPI::testRequest(CommunicationRequest * request) { MPI_Status status; int flag; CommunicationRequestMPI * req_mpi = static_cast(request); MPI_Request * req = req_mpi->getMPIRequest(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Test(req, &flag, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Test."); return (flag != 0); } /* -------------------------------------------------------------------------- */ void StaticCommunicatorMPI::wait(CommunicationRequest * request) { MPI_Status status; CommunicationRequestMPI * req_mpi = static_cast(request); MPI_Request * req = req_mpi->getMPIRequest(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Wait(req, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Wait."); } /* -------------------------------------------------------------------------- */ void StaticCommunicatorMPI::waitAll( std::vector & requests) { MPI_Status status; std::vector::iterator it; for (it = requests.begin(); it != requests.end(); ++it) { MPI_Request * req = static_cast(*it)->getMPIRequest(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Wait(req, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Wait."); } } /* -------------------------------------------------------------------------- */ void StaticCommunicatorMPI::barrier() { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Barrier(communicator); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::reduce(T * values, int nb_values, const SynchronizerOperation & op, int root) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Reduce(MPI_IN_PLACE, values, nb_values, type, MPITypeWrapper::getMPISynchronizerOperation(op), root, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allreduce."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::allReduce(T * values, int nb_values, const SynchronizerOperation & op) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Allreduce(MPI_IN_PLACE, values, nb_values, type, MPITypeWrapper::getMPISynchronizerOperation(op), communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allreduce."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::allGather(T * values, int nb_values) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Allgather(MPI_IN_PLACE, nb_values, type, values, nb_values, type, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allgather."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::allGatherV(T * values, int * nb_values) { MPI_Comm communicator = mpi_data->getMPICommunicator(); int * displs = new int[psize]; displs[0] = 0; for (int i = 1; i < psize; ++i) { displs[i] = displs[i - 1] + nb_values[i - 1]; } MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Allgatherv(MPI_IN_PLACE, *nb_values, type, values, nb_values, displs, type, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Gather."); delete[] displs; } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::gather(T * values, int nb_values, int root) { MPI_Comm communicator = mpi_data->getMPICommunicator(); T * send_buf = NULL, * recv_buf = NULL; if (prank == root) { send_buf = (T *)MPI_IN_PLACE; recv_buf = values; } else { send_buf = values; } MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Gather(send_buf, nb_values, type, recv_buf, nb_values, type, root, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Gather."); } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::gatherV(T * values, int * nb_values, int root) { MPI_Comm communicator = mpi_data->getMPICommunicator(); int * displs = NULL; if (prank == root) { displs = new int[psize]; displs[0] = 0; for (int i = 1; i < psize; ++i) { displs[i] = displs[i - 1] + nb_values[i - 1]; } } T * send_buf = NULL, * recv_buf = NULL; if (prank == root) { send_buf = (T *)MPI_IN_PLACE; recv_buf = values; } else send_buf = values; MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Gatherv(send_buf, *nb_values, type, recv_buf, nb_values, displs, type, root, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Gather."); if (prank == root) { delete[] displs; } } /* -------------------------------------------------------------------------- */ template void StaticCommunicatorMPI::broadcast(T * values, int nb_values, int root) { MPI_Comm communicator = mpi_data->getMPICommunicator(); MPI_Datatype type = MPITypeWrapper::getMPIDatatype(); #if !defined(AKANTU_NDEBUG) int ret = #endif MPI_Bcast(values, nb_values, type, root, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Gather."); } /* -------------------------------------------------------------------------- */ int StaticCommunicatorMPI::getMaxTag() { return MPI_TAG_UB; } /* -------------------------------------------------------------------------- */ int StaticCommunicatorMPI::getMinTag() { return 0; } /* -------------------------------------------------------------------------- */ // template // MPI_Datatype StaticCommunicatorMPI::getMPIDatatype() { // return MPI_DATATYPE_NULL; // } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_CHAR; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_FLOAT; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_DOUBLE; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_LONG_DOUBLE; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_INT; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_UNSIGNED; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_LONG; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_UNSIGNED_LONG; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_LONG_LONG; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype() { return MPI_UNSIGNED_LONG_LONG; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype >() { return MPI_DOUBLE_INT; } template <> MPI_Datatype MPITypeWrapper::getMPIDatatype >() { return MPI_FLOAT_INT; } /* -------------------------------------------------------------------------- */ /* Template instantiation */ /* -------------------------------------------------------------------------- */ #define AKANTU_MPI_COMM_INSTANTIATE(T) \ template void StaticCommunicatorMPI::send(T * buffer, Int size, \ Int receiver, Int tag); \ template void StaticCommunicatorMPI::receive(T * buffer, Int size, \ Int sender, Int tag); \ template CommunicationRequest * StaticCommunicatorMPI::asyncSend( \ T * buffer, Int size, Int receiver, Int tag); \ template CommunicationRequest * StaticCommunicatorMPI::asyncReceive( \ T * buffer, Int size, Int sender, Int tag); \ template void StaticCommunicatorMPI::probe(Int sender, Int tag, \ CommunicationStatus & status); \ template void StaticCommunicatorMPI::allGather(T * values, \ int nb_values); \ template void StaticCommunicatorMPI::allGatherV(T * values, \ int * nb_values); \ template void StaticCommunicatorMPI::gather(T * values, int nb_values, \ int root); \ template void StaticCommunicatorMPI::gatherV(T * values, int * nb_values, \ int root); \ template void StaticCommunicatorMPI::broadcast(T * values, int nb_values, \ int root); \ template void StaticCommunicatorMPI::allReduce( \ T * values, int nb_values, const SynchronizerOperation & op); AKANTU_MPI_COMM_INSTANTIATE(Real); AKANTU_MPI_COMM_INSTANTIATE(UInt); AKANTU_MPI_COMM_INSTANTIATE(Int); AKANTU_MPI_COMM_INSTANTIATE(char); template void StaticCommunicatorMPI::send >( SCMinMaxLoc * buffer, Int size, Int receiver, Int tag); template void StaticCommunicatorMPI::receive >( SCMinMaxLoc * buffer, Int size, Int sender, Int tag); template CommunicationRequest * StaticCommunicatorMPI::asyncSend >( SCMinMaxLoc * buffer, Int size, Int receiver, Int tag); template CommunicationRequest * StaticCommunicatorMPI::asyncReceive >( SCMinMaxLoc * buffer, Int size, Int sender, Int tag); template void StaticCommunicatorMPI::probe >( Int sender, Int tag, CommunicationStatus & status); template void StaticCommunicatorMPI::allGather >( SCMinMaxLoc * values, int nb_values); template void StaticCommunicatorMPI::allGatherV >( SCMinMaxLoc * values, int * nb_values); template void StaticCommunicatorMPI::gather >( SCMinMaxLoc * values, int nb_values, int root); template void StaticCommunicatorMPI::gatherV >( SCMinMaxLoc * values, int * nb_values, int root); template void StaticCommunicatorMPI::broadcast >( SCMinMaxLoc * values, int nb_values, int root); template void StaticCommunicatorMPI::allReduce >( SCMinMaxLoc * values, int nb_values, const SynchronizerOperation & op); #if AKANTU_INTEGER_SIZE > 4 AKANTU_MPI_COMM_INSTANTIATE(int); #endif __END_AKANTU__ diff --git a/src/synchronizer/static_communicator_mpi.hh b/src/synchronizer/static_communicator_mpi.hh index 53c56c833..8186ce8a1 100644 --- a/src/synchronizer/static_communicator_mpi.hh +++ b/src/synchronizer/static_communicator_mpi.hh @@ -1,129 +1,128 @@ /** * @file static_communicator_mpi.hh * * @author Nicolas Richart * * @date creation: Sun Sep 05 2010 * @date last modification: Wed Jan 13 2016 * * @brief class handling parallel communication trough MPI * * @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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STATIC_COMMUNICATOR_MPI_HH__ #define __AKANTU_STATIC_COMMUNICATOR_MPI_HH__ /* -------------------------------------------------------------------------- */ #include "real_static_communicator.hh" /* -------------------------------------------------------------------------- */ #include - __BEGIN_AKANTU__ class MPITypeWrapper; /* -------------------------------------------------------------------------- */ class StaticCommunicatorMPI : public RealStaticCommunicator { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - StaticCommunicatorMPI(int & argc, char ** & argv); + StaticCommunicatorMPI(int & argc, char **& argv); virtual ~StaticCommunicatorMPI(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: + template void send(T * buffer, Int size, Int receiver, Int tag); + template void receive(T * buffer, Int size, Int sender, Int tag); - template void send (T * buffer, Int size, Int receiver, Int tag); - template void receive(T * buffer, Int size, Int sender, Int tag); - - template CommunicationRequest * asyncSend (T * buffer, Int size, - Int receiver, Int tag); - template CommunicationRequest * asyncReceive(T * buffer, Int size, - Int sender, Int tag); + template + CommunicationRequest * asyncSend(T * buffer, Int size, Int receiver, Int tag); + template + CommunicationRequest * asyncReceive(T * buffer, Int size, Int sender, + Int tag); - template void probe(Int sender, Int tag, - CommunicationStatus & status); + template + void probe(Int sender, Int tag, CommunicationStatus & status); - template void allGather (T * values, int nb_values); - template void allGatherV(T * values, int * nb_values); + template void allGather(T * values, int nb_values); + template void allGatherV(T * values, int * nb_values); - template void gather (T * values, int nb_values, int root); - template void gatherV(T * values, int * nb_values, int root); + template void gather(T * values, int nb_values, int root); + template void gatherV(T * values, int * nb_values, int root); - template void broadcast(T * values, int nb_values, int root); + template void broadcast(T * values, int nb_values, int root); bool testRequest(CommunicationRequest * request); - void wait (CommunicationRequest * request); + void wait(CommunicationRequest * request); void waitAll(std::vector & requests); void barrier(); - template void reduce(T * values, int nb_values, - const SynchronizerOperation & op, - int root); - template void allReduce(T * values, int nb_values, - const SynchronizerOperation & op); + template + void reduce(T * values, int nb_values, const SynchronizerOperation & op, + int root); + template + void allReduce(T * values, int nb_values, const SynchronizerOperation & op); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: const MPITypeWrapper & getMPITypeWrapper() const { return *mpi_data; } MPITypeWrapper & getMPITypeWrapper() { return *mpi_data; } int getMinTag(); int getMaxTag(); private: void setRank(int prank) { this->prank = prank; } void setSize(int psize) { this->psize = psize; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: friend class MPITypeWrapper; MPITypeWrapper * mpi_data; + bool is_externaly_initialized; }; - /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "static_communicator_mpi_inline_impl.hh" __END_AKANTU__ #endif /* __AKANTU_STATIC_COMMUNICATOR_MPI_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local/CMakeLists.txt index f807bd4ec..332dcd713 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local/CMakeLists.txt @@ -1,52 +1,52 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # @author Clement Roux # # @date creation: Fri Sep 03 2010 # @date last modification: Mon Dec 07 2015 # # @brief configuration for materials tests # # @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 # #=============================================================================== add_mesh(test_material_non_local_mesh mesh.geo 2 1 OUTPUT mesh.msh) add_mesh(test_material_damage_non_local_mesh mesh_section_gap.geo 2 1 OUTPUT mesh_section_gap.msh) register_test(test_material_damage_non_local SOURCES test_material_damage_non_local.cc DEPENDS test_material_damage_non_local_mesh FILES_TO_COPY material_damage_non_local.dat DIRECTORIES_TO_CREATE paraview PACKAGE damage_non_local ) register_test(test_material_non_local - SOURCES test_material_non_local.cc custom_non_local_test_material.cc + SOURCES test_material_non_local.cc custom_non_local_test_material.cc custom_non_local_test_material.hh DEPENDS test_material_non_local_mesh FILES_TO_COPY material.dat PACKAGE damage_non_local )