diff --git a/cmake/libakantu/v2/printers.py b/cmake/libakantu/v2/printers.py index 1dae8062a..4492e2d67 100755 --- a/cmake/libakantu/v2/printers.py +++ b/cmake/libakantu/v2/printers.py @@ -1,218 +1,179 @@ #!/usr/bin/env python # encoding: utf-8 # -# Inspired from boost's pretty printers from Rüdiger Sonderfeld +# Inspired from boost's pretty printers from +# Rüdiger Sonderfeld # and from Pretty-printers for libstc++ from Free Software Foundation, Inc. # import gdb -import itertools -import numpy import re -#import libstdcxx.v6.printers as std +# import libstdcxx.v6.printers as std __use_gdb_pp__ = True try: import gdb.printing except ImportError: __use_gdb_pp__ = False class AkantuPrinter(object): regex = None @classmethod def supports(cls, typename): return cls.regex.search(typename) @classmethod def get_basic_type(cls, value): """ Determines the type associated to a value""" _type = value.type # If it points to a reference, get the reference. if _type.code == gdb.TYPE_CODE_REF: - _type = _type.target () + _type = _type.target() # Get the unqualified type, stripped of typedefs. _type = _type.unqualified().strip_typedefs() return _type.tag if __use_gdb_pp__: - __akantu_pretty_printers__ = gdb.printing.RegexpCollectionPrettyPrinter("libakantu-v2") + __akantu_pretty_printers__ = \ + gdb.printing.RegexpCollectionPrettyPrinter("libakantu-v2") else: class AkantuPrettyPrinters(object): def __init__(self, name): super(AkantuPrettyPrinters, self).__init__() self.name = name self.printers = {} def add_printer(self, name, regex, printer): self.printers[name] = printer def __call__(self, val): typename = AkantuPrinter.get_basic_type(val) if not typename: return None for name, printer in self.printers.iteritems(): if(printer.supports(typename)): return printer return None __akantu_pretty_printers__ = AkantuPrettyPrinters("libakantu-v2") def register_pretty_printer(pretty_printer): "Registers a Pretty Printer" __akantu_pretty_printers__.add_printer(pretty_printer.name, pretty_printer.regex, pretty_printer) return pretty_printer -@register_pretty_printer + +# @register_pretty_printer class AkaTensorPrinter(AkantuPrinter): """Pretty printer for akantu::Tensor""" regex = re.compile('^akantu::Tensor<(.*), +(.*), +(.*)>$') name = 'akantu::Tensor' value = None typename = "" ptr = None dims = [] ndims = 0 def pretty_print(self): def ij2str(i, j, m): return "{0}".format((self.ptr+m*j + i).dereference()) + def line(i, m, n): - return "[{0}]".format(", ".join((ij2str(i,j, m) for j in range(n)))) + return "[{0}]".format(", ".join((ij2str(i, j, m) for j in + range(n)))) m = int(self.dims[0]) if (self.ndims == 1): n = 1 else: n = int(self.dims[1]) return "[{0}]".format(", ".join(line(i, m, n) for i in range(m))) def __init__(self, value): self.typename = self.get_basic_type(value) self.value = value self.ptr = self.value['values'] self.dims = self.value['n'] def children(self): yield ('values', self.pretty_print()) yield ('wrapped', self.value['wrapped']) @register_pretty_printer class AkaVectorPrinter(AkaTensorPrinter): """Pretty printer for akantu::Vector""" regex = re.compile('^akantu::Vector<(.*)>$') name = 'akantu::Vector' n = 0 ptr = 0x0 def __init__(self, value): super(AkaVectorPrinter, self).__init__(value) self.ndims = 1 - def to_string (self): + def to_string(self): m = self.regex.search(self.typename) return 'Vector<{0}>({1}) [{2}]'.format(m.group(1), int(self.dims[0]), str(self.ptr)) + @register_pretty_printer class AkaMatrixPrinter(AkaTensorPrinter): """Pretty printer for akantu::Matrix""" - regex = re.compile('^akantu::Matrix<(.*)>$'); - name = 'akantu::Matrix'; + regex = re.compile('^akantu::Matrix<(.*)>$') + name = 'akantu::Matrix' def __init__(self, value): super(AkaMatrixPrinter, self).__init__(value) self.ndims = 2 - def to_string (self): + def to_string(self): m = self.regex.search(self.typename) return 'Matrix<%s>(%d, %d) [%s]' % (m.group(1), int(self.dims[0]), int(self.dims[1]), str(self.ptr)) + @register_pretty_printer class AkaElementPrinter(AkantuPrinter): """Pretty printer for akantu::Element""" regex = re.compile('^akantu::Element$') name = 'akantu::Element' def __init__(self, value): self.typename = self.get_basic_type(value) self.value = value - self.element = self.value['element'] - self.eltype = self.value['type'] + self.element = self.value['element'] + self.eltype = self.value['type'] self.ghost_type = self.value['ghost_type'] - def to_string (self): - m = self.regex.search(self.typename) - return 'Element({0}, {1}, {2})'.format(self.element, self.eltype, self.ghost_type) + def to_string(self): + return 'Element({0}, {1}, {2})'.format(self.element, self.eltype, + self.ghost_type) -#@register_pretty_printer -#class AkantuElementTypeMapArrayPrinter(std.StdMapPrinter): -# """Pretty printer for akantu::ElementTypeMapArray""" -# regex = re.compile('^akantu::ElementTypeMapArray<(.*), akantu::(.*)>$') -# name = 'akantu::ElementTypeMapArray' -# -# def __init__ (self, value): -# self.typename = AkantuPrinter.get_basic_type(value) -# self.value = value -# -# def to_string (self): -# m = self.regex.search(self.typename) -# return '{0}MapArray<{1}> with {2} elements and {3} ghost elements'.format(m.group(2), -# m.group(1), -# len (std.RbtreeIterator (self.value['data'])), -# len (std.RbtreeIterator (self.value['data']))) -# -# def children (self): -# yield ('elements:', self.value['data'].children()) -# yield ('ghost elements:', self.value['data'].children()) -# -# def display_hint (self): -# return 'map' -# -#@register_pretty_printer -#class AkantuElementTypeMapPrinter(std.StdMapPrinter): -# """Pretty printer for akantu::ElementTypeMap""" -# regex = re.compile('^akantu::ElementTypeMap<(.*), akantu::(.*)>$') -# name = 'akantu::ElementTypeMap' -# -# def __init__ (self, value): -# self.typename = AkantuPrinter.get_basic_type(value) -# self.value = value -# -# def to_string (self): -# m = self.regex.search(self.typename) -# return '%s<%s> with %d elements' % (m.group(2), m.group(1), -# len (std.RbtreeIterator (self.value['data']))) -# -# def display_hint (self): -# return 'map' -# def register_akantu_printers(obj): "Register Akantu Pretty Printers." if __use_gdb_pp__: gdb.printing.register_pretty_printer(obj, __akantu_pretty_printers__) else: if obj is None: obj = gdb obj.pretty_printers.append(__akantu_pretty_printers__) - - diff --git a/src/common/aka_common_inline_impl.cc b/src/common/aka_common_inline_impl.cc index fda1c4038..33594719f 100644 --- a/src/common/aka_common_inline_impl.cc +++ b/src/common/aka_common_inline_impl.cc @@ -1,433 +1,438 @@ /** * @file aka_common_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Oct 15 2015 * * @brief inline implementations of common akantu type descriptions * * @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 * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /// standard output stream operator for GhostType inline std::ostream & operator<<(std::ostream & stream, GhostType type) { switch (type) { case _not_ghost: stream << "not_ghost"; break; case _ghost: stream << "ghost"; break; case _casper: stream << "Casper the friendly ghost"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for TimeStepSolverType inline std::ostream & operator<<(std::ostream & stream, const TimeStepSolverType & type) { switch (type) { case _tsst_static: stream << "static"; break; case _tsst_dynamic: stream << "dynamic"; break; case _tsst_dynamic_lumped: stream << "dynamic_lumped"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard input stream operator for TimeStepSolverType inline std::istream & operator>>(std::istream & stream, TimeStepSolverType & type) { std::string str; stream >> str; if (str == "static") type = _tsst_static; else if (str == "dynamic") type = _tsst_dynamic; else if (str == "dynamic_lumped") type = _tsst_dynamic_lumped; else { + AKANTU_DEBUG_ERROR("The type " + << str << " is not a recognized TimeStepSolverType"); + stream.setstate(std::ios::failbit); } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for NonLinearSolverType inline std::ostream & operator<<(std::ostream & stream, const NonLinearSolverType & type) { switch (type) { case _nls_linear: stream << "linear"; break; case _nls_newton_raphson: stream << "newton_raphson"; break; case _nls_newton_raphson_modified: stream << "newton_raphson_modified"; break; case _nls_lumped: stream << "lumped"; break; case _nls_auto: stream << "auto"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard input stream operator for NonLinearSolverType inline std::istream & operator>>(std::istream & stream, NonLinearSolverType & type) { std::string str; stream >> str; if (str == "linear") type = _nls_linear; else if (str == "newton_raphson") type = _nls_newton_raphson; else if (str == "newton_raphson_modified") type = _nls_newton_raphson_modified; else if (str == "lumped") type = _nls_lumped; else if (str == "auto") type = _nls_auto; else type = _nls_auto; return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for IntegrationSchemeType inline std::ostream & operator<<(std::ostream & stream, const IntegrationSchemeType & type) { switch (type) { case _ist_pseudo_time: stream << "pseudo_time"; break; case _ist_forward_euler: stream << "forward_euler"; break; case _ist_trapezoidal_rule_1: stream << "trapezoidal_rule_1"; break; case _ist_backward_euler: stream << "backward_euler"; break; case _ist_central_difference: stream << "central_difference"; break; case _ist_fox_goodwin: stream << "fox_goodwin"; break; case _ist_trapezoidal_rule_2: stream << "trapezoidal_rule_2"; break; case _ist_linear_acceleration: stream << "linear_acceleration"; break; case _ist_newmark_beta: stream << "newmark_beta"; break; case _ist_generalized_trapezoidal: stream << "generalized_trapezoidal"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard input stream operator for IntegrationSchemeType inline std::istream & operator>>(std::istream & stream, IntegrationSchemeType & type) { std::string str; stream >> str; if (str == "pseudo_time") type = _ist_pseudo_time; else if (str == "forward_euler") type = _ist_forward_euler; else if (str == "trapezoidal_rule_1") type = _ist_trapezoidal_rule_1; else if (str == "backward_euler") type = _ist_backward_euler; else if (str == "central_difference") type = _ist_central_difference; else if (str == "fox_goodwin") type = _ist_fox_goodwin; else if (str == "trapezoidal_rule_2") type = _ist_trapezoidal_rule_2; else if (str == "linear_acceleration") type = _ist_linear_acceleration; else if (str == "newmark_beta") type = _ist_newmark_beta; else if (str == "generalized_trapezoidal") type = _ist_generalized_trapezoidal; else { + AKANTU_DEBUG_ERROR("The type " + << str << " is not a recognized IntegrationSchemeType"); stream.setstate(std::ios::failbit); } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for SynchronizationTag inline std::ostream & operator<<(std::ostream & stream, SynchronizationTag type) { switch (type) { case _gst_smm_mass: stream << "_gst_smm_mass"; break; case _gst_smm_for_gradu: stream << "_gst_smm_for_gradu"; break; case _gst_smm_boundary: stream << "_gst_smm_boundary"; break; case _gst_smm_uv: stream << "_gst_smm_uv"; break; case _gst_smm_res: stream << "_gst_smm_res"; break; case _gst_smm_init_mat: stream << "_gst_smm_init_mat"; break; case _gst_smm_stress: stream << "_gst_smm_stress"; break; case _gst_smmc_facets: stream << "_gst_smmc_facets"; break; case _gst_smmc_facets_conn: stream << "_gst_smmc_facets_conn"; break; case _gst_smmc_facets_stress: stream << "_gst_smmc_facets_stress"; break; case _gst_smmc_damage: stream << "_gst_smmc_damage"; break; case _gst_giu_global_conn: stream << "_gst_giu_global_conn"; break; case _gst_ce_groups: stream << "_gst_ce_groups"; break; case _gst_gm_clusters: stream << "_gst_gm_clusters"; break; case _gst_htm_capacity: stream << "_gst_htm_capacity"; break; case _gst_htm_temperature: stream << "_gst_htm_temperature"; break; case _gst_htm_gradient_temperature: stream << "_gst_htm_gradient_temperature"; break; case _gst_htm_phi: stream << "_gst_htm_phi"; break; case _gst_htm_gradient_phi: stream << "_gst_htm_gradient_phi"; break; case _gst_mnl_for_average: stream << "_gst_mnl_for_average"; break; case _gst_mnl_weight: stream << "_gst_mnl_weight"; break; case _gst_nh_criterion: stream << "_gst_nh_criterion"; break; case _gst_test: stream << "_gst_test"; break; case _gst_user_1: stream << "_gst_user_1"; break; case _gst_user_2: stream << "_gst_user_2"; break; case _gst_material_id: stream << "_gst_material_id"; break; case _gst_for_dump: stream << "_gst_for_dump"; break; case _gst_cf_nodal: stream << "_gst_cf_nodal"; break; case _gst_cf_incr: stream << "_gst_cf_incr"; break; case _gst_solver_solution: stream << "_gst_solver_solution"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for SolveConvergenceCriteria inline std::ostream & operator<<(std::ostream & stream, const SolveConvergenceCriteria & criteria) { switch (criteria) { case _scc_residual: stream << "_scc_residual"; break; case _scc_solution: stream << "_scc_solution"; break; case _scc_residual_mass_wgh: stream << "_scc_residual_mass_wgh"; break; } return stream; } inline std::istream & operator>>(std::istream & stream, SolveConvergenceCriteria & criteria) { std::string str; stream >> str; if (str == "residual") criteria = _scc_residual; else if (str == "solution") criteria = _scc_solution; else if (str == "residual_mass_wgh") criteria = _scc_residual_mass_wgh; else { stream.setstate(std::ios::failbit); } return stream; } /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str) { std::string lstr = str; std::transform(lstr.begin(), lstr.end(), lstr.begin(), (int (*)(int))tolower); return lstr; } /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim) { std::string trimed = to_trim; // left trim trimed.erase(trimed.begin(), std::find_if(trimed.begin(), trimed.end(), std::not1(std::ptr_fun(isspace)))); // right trim trimed.erase(std::find_if(trimed.rbegin(), trimed.rend(), std::not1(std::ptr_fun(isspace))).base(), trimed.end()); return trimed; } __END_AKANTU__ #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template std::string printMemorySize(UInt size) { Real real_size = size * sizeof(T); UInt mult = (std::log(real_size) / std::log(2)) / 10; std::stringstream sstr; real_size /= Real(1 << (10 * mult)); sstr << std::setprecision(2) << std::fixed << real_size; std::string size_prefix; switch (mult) { case 0: sstr << ""; break; case 1: sstr << "Ki"; break; case 2: sstr << "Mi"; break; case 3: sstr << "Gi"; break; // I started on this type of machines // (32bit computers) (Nicolas) case 4: sstr << "Ti"; break; case 5: sstr << "Pi"; break; case 6: sstr << "Ei"; break; // theoritical limit of RAM of the current // computers in 2014 (64bit computers) (Nicolas) case 7: sstr << "Zi"; break; case 8: sstr << "Yi"; break; default: AKANTU_DEBUG_ERROR( "The programmer in 2014 didn't thought so far (even wikipedia does not " "go further)." << " You have at least 1024 times more than a yobibit of RAM!!!" << " Just add the prefix corresponding in this switch case."); } sstr << "Byte"; return sstr.str(); } __END_AKANTU__ diff --git a/src/common/aka_error.hh b/src/common/aka_error.hh index 03c85e3f4..102509abf 100644 --- a/src/common/aka_error.hh +++ b/src/common/aka_error.hh @@ -1,338 +1,339 @@ /** * @file aka_error.hh * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Jul 13 2015 * * @brief error management and internal exceptions * * @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_ERROR_HH__ #define __AKANTU_ERROR_HH__ /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ enum DebugLevel { dbl0 = 0, dblError = 0, dblAssert = 0, dbl1 = 1, dblException = 1, dblCritical = 1, dbl2 = 2, dblMajor = 2, dbl3 = 3, dblCall = 3, dblSecondary = 3, dblHead = 3, dbl4 = 4, dblWarning = 4, dbl5 = 5, dblInfo = 5, dbl6 = 6, dblIn = 6, dblOut = 6, dbl7 = 7, dbl8 = 8, dblTrace = 8, dbl9 = 9, dblAccessory = 9, dbl10 = 10, dblDebug = 42, dbl100 = 100, dblDump = 100, dblTest = 1337 }; /* -------------------------------------------------------------------------- */ #define AKANTU_LOCATION \ "(" << __func__ << "(): " << __FILE__ << ":" << __LINE__ << ")" /* -------------------------------------------------------------------------- */ namespace debug { void setDebugLevel(const DebugLevel & level); const DebugLevel & getDebugLevel(); void initSignalHandler(); std::string demangle(const char * symbol); std::string exec(std::string cmd); void printBacktrace(int sig); void exit(int status) __attribute__((noreturn)); /* ------------------------------------------------------------------------ */ /// exception class that can be thrown by akantu class Exception : public std::exception { /* ---------------------------------------------------------------------- */ /* Constructors/Destructors */ /* ---------------------------------------------------------------------- */ protected: Exception(std::string info = "") : _info(info), _file(""), _line(0) {} public: //! full constructor Exception(std::string info, std::string file, unsigned int line) : _info(info), _file(file), _line(line) {} //! destructor virtual ~Exception() throw(){}; /* ---------------------------------------------------------------------- */ /* Methods */ /* ---------------------------------------------------------------------- */ public: virtual const char * what() const throw() { return _info.c_str(); } virtual const char * info() const throw() { std::stringstream stream; stream << debug::demangle(typeid(*this).name()) << " : " << _info << " [" << _file << ":" << _line << "]"; return stream.str().c_str(); } public: void setInfo(std::string info) { _info = info; } void setFile(std::string file) { _file = file; } void setLine(unsigned int line) { _line = line; } /* ---------------------------------------------------------------------- */ /* Class Members */ /* ---------------------------------------------------------------------- */ private: /// exception description and additionals std::string _info; /// file it is thrown from std::string _file; /// ligne it is thrown from unsigned int _line; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Exception & _this) { stream << _this.what(); return stream; } /* -------------------------------------------------------------------------- */ class Debugger { public: Debugger(); virtual ~Debugger(); void exit(int status) __attribute__((noreturn)); void 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) __attribute__((noreturn)); /*----------------------------------------------------------------------- */ template void throwCustomException(const Except & ex, const std::string & info, const std::string & file, unsigned int line) const throw(Except) __attribute__((noreturn)); /*----------------------------------------------------------------------- */ template void throwCustomException(const Except & ex, const std::string & file, unsigned int line) const throw(Except) __attribute__((noreturn)); void printMessage(const std::string & prefix, const DebugLevel & level, const std::string & info) const; void setOutStream(std::ostream & out) { cout = &out; } std::ostream & getOutStream() { return *cout; } public: void setParallelContext(int rank, int size); void setDebugLevel(const DebugLevel & level); const DebugLevel & getDebugLevel() const; void setLogFile(const std::string & filename); std::ostream & getOutputStream(); inline bool testLevel(const DebugLevel & level) const { return (this->level >= (level)); } void printBacktrace(bool on_off) { this->print_backtrace = on_off; } bool printBacktrace() { return this->print_backtrace; } private: std::string parallel_context; std::ostream * cout; bool file_open; DebugLevel level; bool print_backtrace; }; extern Debugger debugger; } /* -------------------------------------------------------------------------- */ #define AKANTU_STRINGSTREAM_IN(_str, _sstr) \ ; \ do { \ std::stringstream _dbg_s_info; \ _dbg_s_info << _sstr; \ _str = _dbg_s_info.str(); \ } while (0) /* -------------------------------------------------------------------------- */ #define AKANTU_EXCEPTION(info) AKANTU_EXCEPTION_(info, false) #define AKANTU_SILENT_EXCEPTION(info) AKANTU_EXCEPTION_(info, true) #define AKANTU_EXCEPTION_(info, silent) \ do { \ std::stringstream _dbg_str; \ _dbg_str << info; \ std::stringstream _dbg_loc; \ _dbg_loc << AKANTU_LOCATION; \ ::akantu::debug::debugger.throwException( \ _dbg_str.str(), __FILE__, __LINE__, silent, _dbg_loc.str()); \ } while (0) #define AKANTU_CUSTOM_EXCEPTION_INFO(ex, info) \ do { \ std::d::stringstream _dbg_str; \ _dbg_str << info; \ ::akantu::debug::debugger.throwCustomException(ex, _dbg_str.str(), \ __FILE__, __LINE__); \ } while (0) #define AKANTU_CUSTOM_EXCEPTION(ex) \ do { \ ::akantu::debug::debugger.throwCustomException(ex, __FILE__, __LINE__); \ } while (0) /* -------------------------------------------------------------------------- */ #ifdef AKANTU_NDEBUG #define AKANTU_DEBUG_TEST(level) (false) #define AKANTU_DEBUG_LEVEL_IS_TEST() \ (::akantu::debug::debugger.testLevel(dblTest)) #define AKANTU_DEBUG(level, info) #define AKANTU_DEBUG_(pref, level, info) #define AKANTU_DEBUG_IN() #define AKANTU_DEBUG_OUT() #define AKANTU_DEBUG_INFO(info) #define AKANTU_DEBUG_WARNING(info) #define AKANTU_DEBUG_TRACE(info) #define AKANTU_DEBUG_ASSERT(test, info) #define AKANTU_DEBUG_ERROR(info) AKANTU_EXCEPTION(info) #define AKANTU_DEBUG_TO_IMPLEMENT() ::akantu::debug::debugger.exit(EXIT_FAILURE) /* -------------------------------------------------------------------------- */ #else #define AKANTU_DEBUG(level, info) AKANTU_DEBUG_(" ", level, info) #define AKANTU_DEBUG_(pref, level, info) \ do { \ std::string _dbg_str; \ AKANTU_STRINGSTREAM_IN(_dbg_str, info << " " << AKANTU_LOCATION); \ ::akantu::debug::debugger.printMessage(pref, level, _dbg_str); \ } while (0) #define AKANTU_DEBUG_TEST(level) (::akantu::debug::debugger.testLevel(level)) #define AKANTU_DEBUG_LEVEL_IS_TEST() \ (::akantu::debug::debugger.testLevel(dblTest)) #define AKANTU_DEBUG_IN() \ AKANTU_DEBUG_("==>", ::akantu::dblIn, __func__ << "()") #define AKANTU_DEBUG_OUT() \ AKANTU_DEBUG_("<==", ::akantu::dblOut, __func__ << "()") #define AKANTU_DEBUG_INFO(info) AKANTU_DEBUG_("---", ::akantu::dblInfo, info) #define AKANTU_DEBUG_WARNING(info) \ AKANTU_DEBUG_("/!\\", ::akantu::dblWarning, info) #define AKANTU_DEBUG_TRACE(info) AKANTU_DEBUG_(">>>", ::akantu::dblTrace, info) #define AKANTU_DEBUG_ASSERT(test, info) \ do { \ if (!(test)) { \ AKANTU_DEBUG_("!!! ", ::akantu::dblAssert, "assert [" << #test << "] " \ << info); \ ::akantu::debug::debugger.exit(EXIT_FAILURE); \ } \ } while (0) #define AKANTU_DEBUG_ERROR(info) \ do { \ AKANTU_DEBUG_("!!! ", ::akantu::dblError, info); \ ::akantu::debug::debugger.exit(EXIT_FAILURE); \ } while (0) #define AKANTU_DEBUG_TO_IMPLEMENT() \ AKANTU_DEBUG_ERROR(__func__ << " : not implemented yet !") #endif // AKANTU_NDEBUG /* -------------------------------------------------------------------------- */ namespace debug { /* ---------------------------------------------------------------------- */ template void Debugger::throwCustomException(const Except & ex, const std::string & info, const std::string & file, unsigned int line) const throw(Except) { Except & nc_ex = const_cast(ex); nc_ex.setInfo(info); nc_ex.setFile(file); nc_ex.setLine(line); throw ex; } /* ----------------------------------------------------------------------- */ template void Debugger::throwCustomException(const Except & ex, const std::string & file, unsigned int line) const throw(Except) { Except & nc_ex = const_cast(ex); nc_ex.setFile(file); nc_ex.setLine(line); + throw ex; } } } #endif /* __AKANTU_ERROR_HH__ */ diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh index fb1322bf8..1c4a8ba8a 100644 --- a/src/common/aka_types.hh +++ b/src/common/aka_types.hh @@ -1,1167 +1,1156 @@ /** * @file aka_types.hh * * @author Nicolas Richart * * @date creation: Thu Feb 17 2011 * @date last modification: Fri Jan 22 2016 * * @brief description of the "simple" types * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_error.hh" #include "aka_fwd.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_TYPES_HH__ #define __AKANTU_AKA_TYPES_HH__ __BEGIN_AKANTU__ enum NormType { L_1 = 1, L_2 = 2, L_inf = UInt(-1) }; /** * DimHelper is a class to generalize the setup of a dim array from 3 * values. This gives a common interface in the TensorStorage class * independently of its derived inheritance (Vector, Matrix, Tensor3) * @tparam dim */ template struct DimHelper { static inline void setDims(UInt m, UInt n, UInt p, UInt dims[dim]); }; /* -------------------------------------------------------------------------- */ template <> struct DimHelper<1> { static inline void setDims(UInt m, __attribute__((unused)) UInt n, __attribute__((unused)) UInt p, UInt dims[1]) { dims[0] = m; } }; /* -------------------------------------------------------------------------- */ template <> struct DimHelper<2> { static inline void setDims(UInt m, UInt n, __attribute__((unused)) UInt p, UInt dims[2]) { dims[0] = m; dims[1] = n; } }; /* -------------------------------------------------------------------------- */ template <> struct DimHelper<3> { static inline void setDims(UInt m, UInt n, UInt p, UInt dims[3]) { dims[0] = m; dims[1] = n; dims[2] = p; } }; /* -------------------------------------------------------------------------- */ template class TensorStorage; /* -------------------------------------------------------------------------- */ /* Proxy classes */ /* -------------------------------------------------------------------------- */ /** * @class TensorProxy aka_types.hh - * @desc The TensorProxy class is a proxy class to the TensorStorage it handles the + * @desc The TensorProxy class is a proxy class to the TensorStorage it handles + * the * wrapped case. That is to say if an accessor should give access to a Tensor * wrapped on some data, like the Array::iterator they can return a * TensorProxy that will be automatically transformed as a TensorStorage wrapped * on the same data * @tparam T stored type * @tparam ndim order of the tensor * @tparam RetType real derived type */ template class TensorProxy { protected: TensorProxy(T * data, UInt m, UInt n, UInt p) { DimHelper::setDims(m, n, p, this->n); this->values = data; } TensorProxy(const TensorProxy & other) { this->values = other.storage(); for (UInt i = 0; i < ndim; ++i) this->n[i] = other.n[i]; } inline TensorProxy(const TensorStorage & other); + typedef typename RetType::proxy RetTypeProxy; + public: + operator RetType() { return RetType(static_cast(*this)); } + UInt size(UInt i) const { AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim << " dimensions, not " << (i + 1)); return n[i]; } inline UInt size() const { UInt _size = 1; for (UInt d = 0; d < ndim; ++d) _size *= this->n[d]; return _size; } T * storage() const { return values; } - inline TensorProxy & operator=(const RetType & src) { + inline RetTypeProxy & operator=(const RetType & src) { AKANTU_DEBUG_ASSERT( src.size() == this->size(), "You are trying to copy two tensors with different sizes"); memcpy(this->values, src.storage(), this->size() * sizeof(T)); return *this; } - inline TensorProxy & operator=(const TensorProxy & src) { + inline RetTypeProxy & operator=(const RetTypeProxy & src) { AKANTU_DEBUG_ASSERT( src.size() == this->size(), "You are trying to copy two tensors with different sizes"); memcpy(this->values, src.storage(), this->size() * sizeof(T)); return *this; } + template inline RetTypeProxy & operator*=(const O & o) { + RetType(*this) *= o; + return static_cast(*this); + } + + template inline RetTypeProxy & operator/=(const O & o) { + RetType(*this) /= o; + return static_cast(*this); + } + protected: T * values; UInt n[ndim]; }; /* -------------------------------------------------------------------------- */ template class VectorProxy : public TensorProxy > { typedef TensorProxy > parent; typedef Vector type; public: VectorProxy(T * data, UInt n) : parent(data, n, 0, 0) {} VectorProxy(const VectorProxy & src) : parent(src) {} VectorProxy(const Vector & src) : parent(src) {} - VectorProxy & operator=(const type & src) { - parent::operator=(src); - return *this; - } - VectorProxy & operator=(const VectorProxy & src) { - parent::operator=(src); - return *this; - } - T & operator()(UInt index) { return this->values[index]; }; const T & operator()(UInt index) const { return this->values[index]; }; }; template class MatrixProxy : public TensorProxy > { typedef TensorProxy > parent; typedef Matrix type; public: MatrixProxy(T * data, UInt m, UInt n) : parent(data, m, n, 0) {} MatrixProxy(const MatrixProxy & src) : parent(src) {} MatrixProxy(const type & src) : parent(src) {} - MatrixProxy & operator=(const type & src) { - parent::operator=(src); - return *this; - } - MatrixProxy & operator=(const MatrixProxy & src) { - parent::operator=(src); - return *this; - } }; template class Tensor3Proxy : public TensorProxy > { typedef TensorProxy > parent; typedef Tensor3 type; public: Tensor3Proxy(T * data, UInt m, UInt n, UInt k) : parent(data, m, n, k) {} Tensor3Proxy(const Tensor3Proxy & src) : parent(src) {} Tensor3Proxy(const Tensor3 & src) : parent(src) {} - Tensor3Proxy & operator=(const type & src) { - parent::operator=(src); - return *this; - } - Tensor3Proxy & operator=(const Tensor3Proxy & src) { - parent::operator=(src); - return *this; - } }; /* -------------------------------------------------------------------------- */ /* Tensor base class */ /* -------------------------------------------------------------------------- */ template class TensorStorage { public: typedef T value_type; protected: template void copySize(const TensorType & src) { for (UInt d = 0; d < ndim; ++d) this->n[d] = src.size(d); this->_size = src.size(); } TensorStorage() : values(NULL), wrapped(false) { for (UInt d = 0; d < ndim; ++d) this->n[d] = 0; _size = 0; } TensorStorage(const TensorProxy & proxy) { this->copySize(proxy); this->values = proxy.storage(); this->wrapped = true; } protected: TensorStorage(const TensorStorage & src) {} public: TensorStorage(const TensorStorage & src, bool deep_copy) : values(NULL), wrapped(false) { if (deep_copy) this->deepCopy(src); else this->shallowCopy(src); } protected: TensorStorage(UInt m, UInt n, UInt p, const T & def) { DimHelper::setDims(m, n, p, this->n); this->computeSize(); this->values = new T[this->_size]; this->set(def); this->wrapped = false; } TensorStorage(T * data, UInt m, UInt n, UInt p) { DimHelper::setDims(m, n, p, this->n); this->computeSize(); this->values = data; this->wrapped = true; } public: /* ------------------------------------------------------------------------ */ template inline void shallowCopy(const TensorType & src) { this->copySize(src); if (!this->wrapped) delete[] this->values; this->values = src.storage(); this->wrapped = true; } /* ------------------------------------------------------------------------ */ template inline void deepCopy(const TensorType & src) { this->copySize(src); if (!this->wrapped) delete[] this->values; this->values = new T[this->_size]; memcpy((void *)this->values, (void *)src.storage(), this->_size * sizeof(T)); this->wrapped = false; } virtual ~TensorStorage() { if (!this->wrapped) delete[] this->values; } /* ------------------------------------------------------------------------ */ inline TensorStorage & operator=(const RetType & src) { if (this != &src) { if (this->wrapped) { // this test is not sufficient for Tensor of order higher than 1 AKANTU_DEBUG_ASSERT(this->_size == src.size(), "Tensors of different size"); memcpy((void *)this->values, (void *)src.storage(), this->_size * sizeof(T)); } else { deepCopy(src); } } return *this; } /* ------------------------------------------------------------------------ */ template inline RetType & operator+=(const TensorStorage & other) { T * a = this->storage(); T * b = other.storage(); AKANTU_DEBUG_ASSERT( _size == other.size(), "The two tensors do not have the same size, they cannot be subtracted"); for (UInt i = 0; i < _size; ++i) *(a++) += *(b++); return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ template inline RetType & operator-=(const TensorStorage & other) { T * a = this->storage(); T * b = other.storage(); AKANTU_DEBUG_ASSERT( _size == other.size(), "The two tensors do not have the same size, they cannot be subtracted"); for (UInt i = 0; i < _size; ++i) *(a++) -= *(b++); return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ inline RetType & operator+=(const T & x) { T * a = this->values; for (UInt i = 0; i < _size; ++i) *(a++) += x; return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ inline RetType & operator-=(const T & x) { T * a = this->values; for (UInt i = 0; i < _size; ++i) *(a++) -= x; return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ inline RetType & operator*=(const T & x) { T * a = this->storage(); for (UInt i = 0; i < _size; ++i) *(a++) *= x; return *(static_cast(this)); } /* ---------------------------------------------------------------------- */ inline RetType & operator/=(const T & x) { T * a = this->values; for (UInt i = 0; i < _size; ++i) *(a++) /= x; return *(static_cast(this)); } /// Y = \alpha X + Y inline RetType & aXplusY(const TensorStorage & other, const T & alpha = 1.) { AKANTU_DEBUG_ASSERT( _size == other.size(), "The two tensors do not have the same size, they cannot be subtracted"); Math::aXplusY(this->_size, alpha, other.storage(), this->storage()); return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ T * storage() const { return values; } UInt size() const { return _size; } UInt size(UInt i) const { AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim << " dimensions, not " << (i + 1)); return n[i]; }; /* ------------------------------------------------------------------------ */ inline void clear() { memset(values, 0, _size * sizeof(T)); }; inline void set(const T & t) { std::fill_n(values, _size, t); }; template inline void copy(const TensorType & other) { AKANTU_DEBUG_ASSERT( _size == other.size(), "The two tensors do not have the same size, they cannot be copied"); memcpy(values, other.storage(), _size * sizeof(T)); } bool isWrapped() const { return this->wrapped; } protected: friend class Array; inline void computeSize() { _size = 1; for (UInt d = 0; d < ndim; ++d) _size *= this->n[d]; } protected: template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm += std::pow(std::abs(*it), norm_type); return std::pow(_norm, 1. / norm_type); } }; template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm += std::abs(*it); return _norm; } }; template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm += *it * *it; return sqrt(_norm); } }; template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm = std::max(std::abs(*it), _norm); return _norm; } }; public: /*----------------------------------------------------------------------- */ /// "Entrywise" norm norm @f[ \|\boldsymbol{T}\|_p = \left( /// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}} /// @f] template inline T norm() const { return NormHelper::norm(*this); } protected: UInt n[ndim]; UInt _size; T * values; bool wrapped; }; template inline TensorProxy::TensorProxy( const TensorStorage & other) { this->values = other.storage(); for (UInt i = 0; i < ndim; ++i) this->n[i] = other.size(i); } /* -------------------------------------------------------------------------- */ /* Vector */ /* -------------------------------------------------------------------------- */ template class Vector : public TensorStorage > { typedef TensorStorage > parent; public: typedef typename parent::value_type value_type; typedef VectorProxy proxy; public: Vector() : parent() {} Vector(UInt n, const T & def = T()) : parent(n, 0, 0, def) {} Vector(T * data, UInt n) : parent(data, n, 0, 0) {} Vector(const Vector & src, bool deep_copy = true) : parent(src, deep_copy) {} Vector(const VectorProxy & src) : parent(src) {} public: virtual ~Vector(){}; /* ------------------------------------------------------------------------ */ inline Vector & operator=(const Vector & src) { parent::operator=(src); return *this; } /* ------------------------------------------------------------------------ */ inline T & operator()(UInt i) { AKANTU_DEBUG_ASSERT((i < this->n[0]), "Access out of the vector! " << "Index (" << i << ") is out of the vector of size (" << this->n[0] << ")"); return *(this->values + i); } inline const T & operator()(UInt i) const { AKANTU_DEBUG_ASSERT((i < this->n[0]), "Access out of the vector! " << "Index (" << i << ") is out of the vector of size (" << this->n[0] << ")"); return *(this->values + i); } inline T & operator[](UInt i) { return this->operator()(i); } inline const T & operator[](UInt i) const { return this->operator()(i); } /* ------------------------------------------------------------------------ */ inline Vector & operator*=(Real x) { return parent::operator*=(x); } inline Vector & operator/=(Real x) { return parent::operator/=(x); } - /* ------------------------------------------------------------------------ */ inline Vector & operator*=(const Vector & vect) { T * a = this->storage(); T * b = vect.storage(); for (UInt i = 0; i < this->_size; ++i) *(a++) *= *(b++); return *this; } /* ------------------------------------------------------------------------ */ inline Real dot(const Vector & vect) const { return Math::vectorDot(this->values, vect.storage(), this->_size); } /* ------------------------------------------------------------------------ */ inline Real mean() const { Real mean = 0; T * a = this->storage(); for (UInt i = 0; i < this->_size; ++i) mean += *(a++); return mean / this->_size; } /* ------------------------------------------------------------------------ */ inline Vector & crossProduct(const Vector & v1, const Vector & v2) { AKANTU_DEBUG_ASSERT(this->size() == 3, "crossProduct is only defined in 3D (n=" << this->size() << ")"); AKANTU_DEBUG_ASSERT( this->size() == v1.size() && this->size() == v2.size(), "crossProduct is not a valid operation non matching size vectors"); Math::vectorProduct3(v1.storage(), v2.storage(), this->values); return *this; } /* ------------------------------------------------------------------------ */ inline void solve(Matrix & A, const Vector & b) { AKANTU_DEBUG_ASSERT( this->size() == A.rows() && this->_size = A.cols(), "The size of the solution vector mismatches the size of the matrix"); AKANTU_DEBUG_ASSERT( this->_size == b._size, "The rhs vector has a mismatch in size with the matrix"); Math::solve(this->_size, A.storage(), this->values, b.storage()); } /* ------------------------------------------------------------------------ */ template inline void mul(const Matrix & A, const Vector & x, Real alpha = 1.0); /* ------------------------------------------------------------------------ */ inline Real norm() const { return parent::template norm(); } template inline Real norm() const { return parent::template norm(); } /* ------------------------------------------------------------------------ */ inline void normalize() { Real n = norm(); operator/=(n); } /* ------------------------------------------------------------------------ */ /// norm of (*this - x) inline Real distance(const Vector & y) const { Real * vx = this->values; Real * vy = y.storage(); Real sum_2 = 0; for (UInt i = 0; i < this->_size; ++i, ++vx, ++vy) sum_2 += (*vx - *vy) * (*vx - *vy); return sqrt(sum_2); } /* ------------------------------------------------------------------------ */ inline bool equal(const Vector & v, Real tolerance = Math::getTolerance()) const { T * a = this->storage(); T * b = v.storage(); UInt i = 0; while (i < this->_size && (std::abs(*(a++) - *(b++)) < tolerance)) ++i; return i == this->_size; } /* ------------------------------------------------------------------------ */ inline short compare(const Vector & v, Real tolerance = Math::getTolerance()) const { T * a = this->storage(); T * b = v.storage(); for (UInt i(0); i < this->_size; ++i, ++a, ++b) { if (std::abs(*a - *b) > tolerance) return (((*a - *b) > tolerance) ? 1 : -1); } return 0; } /* ------------------------------------------------------------------------ */ inline bool operator==(const Vector & v) const { return equal(v); } inline bool operator<(const Vector & v) const { return compare(v) == -1; } inline bool operator>(const Vector & v) const { return compare(v) == 1; } /* ------------------------------------------------------------------------ */ /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << "["; for (UInt i = 0; i < this->_size; ++i) { if (i != 0) stream << ", "; stream << this->values[i]; } stream << "]"; } friend class ::akantu::Array; }; typedef Vector RVector; /* ------------------------------------------------------------------------ */ template <> inline bool Vector::equal(const Vector & v, __attribute__((unused)) Real tolerance) const { UInt * a = this->storage(); UInt * b = v.storage(); UInt i = 0; while (i < this->_size && (*(a++) == *(b++))) ++i; return i == this->_size; } /* ------------------------------------------------------------------------ */ /* Matrix */ /* ------------------------------------------------------------------------ */ template class Matrix : public TensorStorage > { typedef TensorStorage > parent; public: typedef typename parent::value_type value_type; typedef MatrixProxy proxy; public: Matrix() : parent() {} Matrix(UInt m, UInt n, const T & def = T()) : parent(m, n, 0, def) {} Matrix(T * data, UInt m, UInt n) : parent(data, m, n, 0) {} Matrix(const Matrix & src, bool deep_copy = true) : parent(src, deep_copy) {} Matrix(const MatrixProxy & src) : parent(src) {} virtual ~Matrix() {} /* ------------------------------------------------------------------------ */ inline Matrix & operator=(const Matrix & src) { parent::operator=(src); return *this; } public: /* ---------------------------------------------------------------------- */ UInt rows() const { return this->n[0]; } UInt cols() const { return this->n[1]; } /* ---------------------------------------------------------------------- */ inline T & at(UInt i, UInt j) { AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])), "Access out of the matrix! " << "Index (" << i << ", " << j << ") is out of the matrix of size (" << this->n[0] << ", " << this->n[1] << ")"); return *(this->values + i + j * this->n[0]); } inline const T & at(UInt i, UInt j) const { AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])), "Access out of the matrix! " << "Index (" << i << ", " << j << ") is out of the matrix of size (" << this->n[0] << ", " << this->n[1] << ")"); return *(this->values + i + j * this->n[0]); } /* ---------------------------------------------------------------------- */ inline T & operator()(UInt i, UInt j) { return this->at(i, j); } inline const T & operator()(UInt i, UInt j) const { return this->at(i, j); } /// give a line vector wrapped on the column i inline VectorProxy operator()(UInt j) { AKANTU_DEBUG_ASSERT(j < this->n[1], "Access out of the matrix! " << "You are trying to access the column vector " << j << " in a matrix of size (" << this->n[0] << ", " << this->n[1] << ")"); return VectorProxy(this->values + j * this->n[0], this->n[0]); } inline const VectorProxy operator()(UInt j) const { AKANTU_DEBUG_ASSERT(j < this->n[1], "Access out of the matrix! " << "You are trying to access the column vector " << j << " in a matrix of size (" << this->n[0] << ", " << this->n[1] << ")"); return VectorProxy(this->values + j * this->n[0], this->n[0]); } inline T & operator[](UInt idx) { return *(this->values + idx); }; inline const T & operator[](UInt idx) const { return *(this->values + idx); }; /* ---------------------------------------------------------------------- */ inline Matrix operator*(const Matrix & B) { Matrix C(this->rows(), B.cols()); C.mul(*this, B); return C; } /* ----------------------------------------------------------------------- */ inline Matrix & operator*=(const T & x) { return parent::operator*=(x); } inline Matrix & operator*=(const Matrix & B) { Matrix C(*this); this->mul(C, B); return *this; } /* ---------------------------------------------------------------------- */ template inline void mul(const Matrix & A, const Matrix & B, T alpha = 1.0) { UInt k = A.cols(); if (tr_A) k = A.rows(); #ifndef AKANTU_NDEBUG if (tr_B) { AKANTU_DEBUG_ASSERT(k == B.cols(), "matrices to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->cols() == B.rows(), "matrices to multiply have no fit dimensions"); } else { AKANTU_DEBUG_ASSERT(k == B.rows(), "matrices to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->cols() == B.cols(), "matrices to multiply have no fit dimensions"); } if (tr_A) { AKANTU_DEBUG_ASSERT(this->rows() == A.cols(), "matrices to multiply have no fit dimensions"); } else { AKANTU_DEBUG_ASSERT(this->rows() == A.rows(), "matrices to multiply have no fit dimensions"); } #endif // AKANTU_NDEBUG Math::matMul(this->rows(), this->cols(), k, alpha, A.storage(), B.storage(), 0., this->storage()); } /* ---------------------------------------------------------------------- */ inline void outerProduct(const Vector & A, const Vector & B) { AKANTU_DEBUG_ASSERT( A.size() == this->rows() && B.size() == this->cols(), "A and B are not compatible with the size of the matrix"); for (UInt i = 0; i < this->rows(); ++i) { for (UInt j = 0; j < this->cols(); ++j) { this->values[i + j * this->rows()] += A[i] * B[j]; } } } private: class EigenSorter { public: EigenSorter(const Vector & eigs) : eigs(eigs) {} bool operator()(const UInt & a, const UInt & b) const { return (eigs(a) > eigs(b)); } private: const Vector & eigs; }; public: /* ---------------------------------------------------------------------- */ inline void eig(Vector & eigenvalues, Matrix & eigenvectors) const { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "eig is not a valid operation on a rectangular matrix"); AKANTU_DEBUG_ASSERT(eigenvalues.size() == this->cols(), "eigenvalues should be of size " << this->cols() << "."); #ifndef AKANTU_NDEBUG if (eigenvectors.storage() != NULL) AKANTU_DEBUG_ASSERT((eigenvectors.rows() == eigenvectors.cols()) && (eigenvectors.rows() == this->cols()), "Eigenvectors needs to be a square matrix of size " << this->cols() << " x " << this->cols() << "."); #endif Matrix tmp = *this; Vector tmp_eigs(eigenvalues.size()); Matrix tmp_eig_vects(eigenvectors.rows(), eigenvectors.cols()); if (tmp_eig_vects.rows() == 0 || tmp_eig_vects.cols() == 0) Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage()); else Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage(), tmp_eig_vects.storage()); Vector perm(eigenvalues.size()); for (UInt i = 0; i < perm.size(); ++i) perm(i) = i; std::sort(perm.storage(), perm.storage() + perm.size(), EigenSorter(tmp_eigs)); for (UInt i = 0; i < perm.size(); ++i) eigenvalues(i) = tmp_eigs(perm(i)); if (tmp_eig_vects.rows() != 0 && tmp_eig_vects.cols() != 0) for (UInt i = 0; i < perm.size(); ++i) { for (UInt j = 0; j < eigenvectors.rows(); ++j) { eigenvectors(j, i) = tmp_eig_vects(j, perm(i)); } } } /* ---------------------------------------------------------------------- */ inline void eig(Vector & eigenvalues) const { Matrix empty; eig(eigenvalues, empty); } /* ---------------------------------------------------------------------- */ inline void eye(T alpha = 1.) { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "eye is not a valid operation on a rectangular matrix"); this->clear(); for (UInt i = 0; i < this->cols(); ++i) { this->values[i + i * this->rows()] = alpha; } } /* ---------------------------------------------------------------------- */ static inline Matrix eye(UInt m, T alpha = 1.) { Matrix tmp(m, m); tmp.eye(alpha); return tmp; } /* ---------------------------------------------------------------------- */ inline T trace() const { AKANTU_DEBUG_ASSERT( this->cols() == this->rows(), "trace is not a valid operation on a rectangular matrix"); T trace = 0.; for (UInt i = 0; i < this->rows(); ++i) { trace += this->values[i + i * this->rows()]; } return trace; } /* ---------------------------------------------------------------------- */ inline Matrix transpose() const { Matrix tmp(this->cols(), this->rows()); for (UInt i = 0; i < this->rows(); ++i) { for (UInt j = 0; j < this->cols(); ++j) { tmp(j, i) = operator()(i, j); } } return tmp; } /* ---------------------------------------------------------------------- */ inline void inverse(const Matrix & A) { AKANTU_DEBUG_ASSERT(A.cols() == A.rows(), "inv is not a valid operation on a rectangular matrix"); AKANTU_DEBUG_ASSERT(this->cols() == A.cols(), "the matrix should have the same size as its inverse"); if (this->cols() == 1) *this->values = 1. / *A.storage(); else if (this->cols() == 2) Math::inv2(A.storage(), this->values); else if (this->cols() == 3) Math::inv3(A.storage(), this->values); else Math::inv(this->cols(), A.storage(), this->values); } /* --------------------------------------------------------------------- */ inline T det() const { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "inv is not a valid operation on a rectangular matrix"); if (this->cols() == 1) return *(this->values); else if (this->cols() == 2) return Math::det2(this->values); else if (this->cols() == 3) return Math::det3(this->values); else return Math::det(this->cols(), this->values); } /* --------------------------------------------------------------------- */ inline T doubleDot(const Matrix & other) const { AKANTU_DEBUG_ASSERT( this->cols() == this->rows(), "doubleDot is not a valid operation on a rectangular matrix"); if (this->cols() == 1) return *(this->values) * *(other.storage()); else if (this->cols() == 2) return Math::matrixDoubleDot22(this->values, other.storage()); else if (this->cols() == 3) return Math::matrixDoubleDot33(this->values, other.storage()); else AKANTU_DEBUG_ERROR("doubleDot is not defined for other spatial dimensions" << " than 1, 2 or 3."); return T(); } /* ---------------------------------------------------------------------- */ /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << "["; for (UInt i = 0; i < this->n[0]; ++i) { if (i != 0) stream << ", "; stream << "["; for (UInt j = 0; j < this->n[1]; ++j) { if (j != 0) stream << ", "; stream << operator()(i, j); } stream << "]"; } stream << "]"; }; }; /* ------------------------------------------------------------------------ */ template template inline void Vector::mul(const Matrix & A, const Vector & x, Real alpha) { #ifndef AKANTU_NDEBUG UInt n = x.size(); if (tr_A) { AKANTU_DEBUG_ASSERT(n == A.rows(), "matrix and vector to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->size() == A.cols(), "matrix and vector to multiply have no fit dimensions"); } else { AKANTU_DEBUG_ASSERT(n == A.cols(), "matrix and vector to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->size() == A.rows(), "matrix and vector to multiply have no fit dimensions"); } #endif Math::matVectMul(A.rows(), A.cols(), alpha, A.storage(), x.storage(), 0., this->storage()); } /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Matrix & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Vector & _this) { _this.printself(stream); return stream; } /* ------------------------------------------------------------------------ */ /* Tensor3 */ /* ------------------------------------------------------------------------ */ template class Tensor3 : public TensorStorage > { typedef TensorStorage > parent; public: typedef typename parent::value_type value_type; typedef Tensor3Proxy proxy; public: Tensor3() : parent(){}; Tensor3(UInt m, UInt n, UInt p, const T & def = T()) : parent(m, n, p, def) {} Tensor3(T * data, UInt m, UInt n, UInt p) : parent(data, m, n, p) {} Tensor3(const Tensor3 & src, bool deep_copy = true) : parent(src, deep_copy) {} public: /* ------------------------------------------------------------------------ */ inline Tensor3 & operator=(const Tensor3 & src) { parent::operator=(src); return *this; } /* ---------------------------------------------------------------------- */ inline T & operator()(UInt i, UInt j, UInt k) { AKANTU_DEBUG_ASSERT( (i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]), "Access out of the tensor3! " << "You are trying to access the element " << "(" << i << ", " << j << ", " << k << ") in a tensor of size (" << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")"); return *(this->values + (k * this->n[0] + i) * this->n[1] + j); } inline const T & operator()(UInt i, UInt j, UInt k) const { AKANTU_DEBUG_ASSERT( (i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]), "Access out of the tensor3! " << "You are trying to access the element " << "(" << i << ", " << j << ", " << k << ") in a tensor of size (" << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")"); return *(this->values + (k * this->n[0] + i) * this->n[1] + j); } inline MatrixProxy operator()(UInt k) { AKANTU_DEBUG_ASSERT((k < this->n[2]), "Access out of the tensor3! " << "You are trying to access the slice " << k << " in a tensor3 of size (" << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")"); return MatrixProxy(this->values + k * this->n[0] * this->n[1], this->n[0], this->n[1]); } inline const MatrixProxy operator()(UInt k) const { AKANTU_DEBUG_ASSERT((k < this->n[2]), "Access out of the tensor3! " << "You are trying to access the slice " << k << " in a tensor3 of size (" << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")"); return MatrixProxy(this->values + k * this->n[0] * this->n[1], this->n[0], this->n[1]); } inline MatrixProxy operator[](UInt k) { return MatrixProxy(this->values + k * this->n[0] * this->n[1], this->n[0], this->n[1]); } inline const MatrixProxy operator[](UInt k) const { return MatrixProxy(this->values + k * this->n[0] * this->n[1], this->n[0], this->n[1]); } }; /* -------------------------------------------------------------------------- */ // support operations for the creation of other vectors /* -------------------------------------------------------------------------- */ template Vector operator*(const T & scalar, const Vector & a) { Vector r(a); r *= scalar; return r; } template Vector operator*(const Vector & a, const T & scalar) { Vector r(a); r *= scalar; return r; } template Vector operator/(const Vector & a, const T & scalar) { Vector r(a); r /= scalar; return r; } template Vector operator+(const Vector & a, const Vector & b) { Vector r(a); r += b; return r; } template Vector operator-(const Vector & a, const Vector & b) { Vector r(a); r -= b; return r; } /* -------------------------------------------------------------------------- */ template Matrix operator*(const T & scalar, const Matrix & a) { Matrix r(a); r *= scalar; return r; } template Matrix operator*(const Matrix & a, const T & scalar) { Matrix r(a); r *= scalar; return r; } template Matrix operator/(const Matrix & a, const T & scalar) { Matrix r(a); r /= scalar; return r; } template Matrix operator+(const Matrix & a, const Matrix & b) { Matrix r(a); r += b; return r; } template Matrix operator-(const Matrix & a, const Matrix & b) { Matrix r(a); r -= b; return r; } __END_AKANTU__ #endif /* __AKANTU_AKA_TYPES_HH__ */ diff --git a/src/fe_engine/fe_engine_template_tmpl.hh b/src/fe_engine/fe_engine_template_tmpl.hh index 0671ffcdc..5bc63a4c4 100644 --- a/src/fe_engine/fe_engine_template_tmpl.hh +++ b/src/fe_engine/fe_engine_template_tmpl.hh @@ -1,1636 +1,1641 @@ /** * @file fe_engine_template_tmpl.hh * * @author Guillaume Anciaux * @author Mauro Corrado * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Feb 15 2011 * @date last modification: Thu Nov 19 2015 * * @brief Template implementation of FEEngineTemplate * * @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 "dof_manager.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template