diff --git a/cmake/libakantu/v3/printers.py b/cmake/libakantu/v3/printers.py index eca4d1e2b..93344b834 100755 --- a/cmake/libakantu/v3/printers.py +++ b/cmake/libakantu/v3/printers.py @@ -1,405 +1,431 @@ #!/usr/bin/env python3 # encoding: utf-8 """ printers.py: gdb pretty printers""" __author__ = "Nicolas Richart" __credits__ = [ "Nicolas Richart ", ] -__copyright__ = "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale" \ - " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ - " en Mécanique des Solides)" +__copyright__ = ( + "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale" + " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" + " en Mécanique des Solides)" +) __license__ = "LGPLv3" # # Inspired from boost's pretty printers from # Rüdiger Sonderfeld # and from Pretty-printers for libstc++ from Free Software Foundation, Inc. # import gdb import re import sys import itertools __use_gdb_pp__ = True try: import gdb.printing except ImportError: __use_gdb_pp__ = False class AkantuPrinter(object): regex = None @classmethod def supports(cls, typename): # print('{0} ~= {1}'.format(typename, cls.regex), file=sys.stderr) return cls.regex.search(typename) @classmethod def get_basic_type(cls, value): - """ Determines the type associated to a 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() # 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-v3") + __akantu_pretty_printers__ = gdb.printing.RegexpCollectionPrettyPrinter( + "libakantu-v3" + ) else: + class AkantuPrettyPrinters(object): def __init__(self, name): super(AkantuPrettyPrinters, self).__init__() self.name = name self.printers = {} @property def enabled(self): return True 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.items(): - if(printer.supports(typename)): + if printer.supports(typename): return printer return None __akantu_pretty_printers__ = AkantuPrettyPrinters("libakantu-v3") def register_pretty_printer(pretty_printer): "Registers a Pretty Printer" - __akantu_pretty_printers__.add_printer(pretty_printer.name, - pretty_printer.regex, - pretty_printer) + __akantu_pretty_printers__.add_printer( + pretty_printer.name, pretty_printer.regex, pretty_printer + ) return pretty_printer @register_pretty_printer class AkaArrayPrinter(AkantuPrinter): """Pretty printer for akantu::Array""" - regex = re.compile('^akantu::Array<(.*?), (true|false)>$') - name = 'akantu::Array' + + regex = re.compile("^akantu::Array<(.*?), (true|false)>$") + name = "akantu::Array" def __init__(self, value): self.typename = self.get_basic_type(value) self.value = value - self.ptr = self.value['values'] - self.size = int(self.value['size_']) - self.nb_component = int(self.value['nb_component']) + self.ptr = self.value["values"] + self.size = int(self.value["size_"]) + self.nb_component = int(self.value["nb_component"]) def display_hint(self): - return 'array' + return "array" def to_string(self): m = self.regex.search(self.typename) - return 'Array<{0}>({1}, {2}) stored at {3}'.format( - m.group(1), self.size, self.nb_component, self.ptr) + return "Array<{0}>({1}, {2}) stored at {3}".format( + m.group(1), self.size, self.nb_component, self.ptr + ) def children(self): _ptr = self.ptr for i in range(self.size): - _values = ["{0}".format((_ptr + j).dereference()) - for j in range(self.nb_component)] + _values = [ + "{0}".format((_ptr + j).dereference()) for j in range(self.nb_component) + ] _ptr = _ptr + self.nb_component - yield ('[{0}]'.format(i), - ('{0}' if self.nb_component == 1 else '[{0}]').format( - ', '.join(_values))) + yield ( + "[{0}]".format(i), + ("{0}" if self.nb_component == 1 else "[{0}]").format( + ", ".join(_values) + ), + ) if sys.version_info[0] > 2: # Python 3 stuff Iterator = object else: # Python 2 stuff class Iterator: """Compatibility mixin for iterators Instead of writing next() methods for iterators, write __next__() methods and use this mixin to make them work in Python 2 as well as Python 3. Idea stolen from the "six" documentation: """ + def next(self): return self.__next__() class RbtreeIterator(Iterator): """ Turn an RB-tree-based container (std::map, std::set etc.) into a Python iterable object. """ def __init__(self, rbtree): - self.size = rbtree['_M_t']['_M_impl']['_M_node_count'] - self.node = rbtree['_M_t']['_M_impl']['_M_header']['_M_left'] + self.size = rbtree["_M_t"]["_M_impl"]["_M_node_count"] + self.node = rbtree["_M_t"]["_M_impl"]["_M_header"]["_M_left"] self.count = 0 def __iter__(self): return self def __len__(self): return int(self.size) def __next__(self): if self.count == self.size: raise StopIteration result = self.node self.count = self.count + 1 if self.count < self.size: # Compute the next node. node = self.node - if node.dereference()['_M_right']: - node = node.dereference()['_M_right'] - while node.dereference()['_M_left']: - node = node.dereference()['_M_left'] + if node.dereference()["_M_right"]: + node = node.dereference()["_M_right"] + while node.dereference()["_M_left"]: + node = node.dereference()["_M_left"] else: - parent = node.dereference()['_M_parent'] - while node == parent.dereference()['_M_right']: + parent = node.dereference()["_M_parent"] + while node == parent.dereference()["_M_right"]: node = parent - parent = parent.dereference()['_M_parent'] - if node.dereference()['_M_right'] != parent: + parent = parent.dereference()["_M_parent"] + if node.dereference()["_M_right"] != parent: node = parent self.node = node return result def get_value_from_aligned_membuf(buf, valtype): """Returns the value held in a __gnu_cxx::__aligned_membuf.""" - return buf['_M_storage'].address.cast(valtype.pointer()).dereference() + return buf["_M_storage"].address.cast(valtype.pointer()).dereference() def get_value_from_Rb_tree_node(node): """Returns the value held in an _Rb_tree_node<_Val>""" try: member = node.type.fields()[1].name - if member == '_M_value_field': + if member == "_M_value_field": # C++03 implementation, node contains the value as a member - return node['_M_value_field'] - elif member == '_M_storage': + return node["_M_value_field"] + elif member == "_M_storage": # C++11 implementation, node stores value in __aligned_membuf valtype = node.type.template_argument(0) - return get_value_from_aligned_membuf(node['_M_storage'], valtype) + return get_value_from_aligned_membuf(node["_M_storage"], valtype) except: # noqa: E722 pass raise ValueError("Unsupported implementation for %s" % str(node.type)) def find_type(orig, name): typ = orig.strip_typedefs() while True: # Strip cv-qualifiers. PR 67440. - search = '%s::%s' % (typ.unqualified(), name) + search = "%s::%s" % (typ.unqualified(), name) try: return gdb.lookup_type(search) except RuntimeError: pass # The type was not found, so try the superclass. We only need # to check the first superclass, so we don't bother with # anything fancier here. field = typ.fields()[0] if not field.is_base_class: raise ValueError("Cannot find type %s::%s" % (str(orig), name)) typ = field.type @register_pretty_printer class AkaElementTypeMapArrayPrinter(AkantuPrinter): """Pretty printer for akantu::ElementTypeMap>""" - regex = re.compile('^akantu::ElementTypeMap\*, akantu::(.*?)>$') # noqa: E501,W605 - name = 'akantu::ElementTypeMapArray' + + regex = re.compile( + "^akantu::ElementTypeMap\*, akantu::(.*?)>$" + ) # noqa: E501,W605 + name = "akantu::ElementTypeMapArray" # Turn an RbtreeIterator into a pretty-print iterator. class _rb_iter(Iterator): def __init__(self, rbiter, type, ghost_type): self.rbiter = rbiter self.count = 0 self.type = type self.ghost_type = ghost_type def __iter__(self): return self def __next__(self): if self.count % 2 == 0: n = next(self.rbiter) n = n.cast(self.type).dereference() n = get_value_from_Rb_tree_node(n) self.pair = n - item = "{0}:{1}".format(n['first'], self.ghost_type) + item = "{0}:{1}".format(n["first"], self.ghost_type) else: - item = self.pair['second'].dereference() - result = ('[{0}]'.format(self.count), item) + item = self.pair["second"].dereference() + result = ("[{0}]".format(self.count), item) self.count = self.count + 1 return result def _iter(self, not_ghost, ghost, type): iter_size = (len(not_ghost), len(ghost)) - it = self._rb_iter(not_ghost, type, '_not_ghost') + it = self._rb_iter(not_ghost, type, "_not_ghost") for _ in range(iter_size[0] * 2): yield next(it) - it = self._rb_iter(ghost, type, '_ghost') + it = self._rb_iter(ghost, type, "_ghost") for _ in range(iter_size[1] * 2): yield next(it) raise StopIteration def __init__(self, value): self.typename = self.get_basic_type(value) self.value = value - self.data = self.value['data'] - self.ghost_data = self.value['ghost_data'] + self.data = self.value["data"] + self.ghost_data = self.value["ghost_data"] def to_string(self): m = self.regex.search(self.typename) - return 'ElementTypMapArray<{0}> with '\ - '{1} _not_ghost and {2} _ghost'.format( - m.group(1), len(RbtreeIterator(self.data)), - len(RbtreeIterator(self.ghost_data))) + return "ElementTypMapArray<{0}> with " "{1} _not_ghost and {2} _ghost".format( + m.group(1), + len(RbtreeIterator(self.data)), + len(RbtreeIterator(self.ghost_data)), + ) def children(self): # m = self.regex.search(self.typename) - rep_type = find_type(self.data.type, '_Rep_type') - node = find_type(rep_type, '_Link_type') + rep_type = find_type(self.data.type, "_Rep_type") + node = find_type(rep_type, "_Link_type") node = node.strip_typedefs() return itertools.chain( self._rb_iter(RbtreeIterator(self.data), node, "_not_ghost"), - self._rb_iter(RbtreeIterator(self.ghost_data), node, "_ghost")) + self._rb_iter(RbtreeIterator(self.ghost_data), node, "_ghost"), + ) def display_hint(self): - return 'map' + return "map" # @register_pretty_printer class AkaTensorPrinter(AkantuPrinter): """Pretty printer for akantu::Tensor""" - regex = re.compile('^akantu::Tensor<(.*), +(.*), +(.*)>$') - name = '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()) + 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): + 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'] + self.ptr = self.value["values"] + self.dims = self.value["n"] def children(self): - yield ('values', self.pretty_print()) - yield ('wrapped', self.value['wrapped']) + 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' + + 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): m = self.regex.search(self.typename) - return 'Vector<{0}>({1}) [{2}]'.format(m.group(1), int(self.dims[0]), - str(self.ptr)) + 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): 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)) + 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' + + 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.ghost_type = self.value['ghost_type'] - self._ek_not_defined = gdb.parse_and_eval('akantu::_not_defined') - self._casper = gdb.parse_and_eval('akantu::_casper') - self._max_uint = gdb.parse_and_eval('(akantu::UInt) -1') + self.element = self.value["element"] + self.eltype = self.value["type"] + self.ghost_type = self.value["ghost_type"] + self._ek_not_defined = gdb.parse_and_eval("akantu::_not_defined") + self._casper = gdb.parse_and_eval("akantu::_casper") + self._max_uint = gdb.parse_and_eval("(akantu::Int) -1") def to_string(self): - if (self.element == self._max_uint) and \ - (self.eltype == self._ek_not_defined) and \ - (self.ghost_type == self._casper): - return 'ElementNull' - - return 'Element({0}, {1}, {2})'.format(self.element, self.eltype, - self.ghost_type) + if ( + (self.element == self._max_uint) + and (self.eltype == self._ek_not_defined) + and (self.ghost_type == self._casper) + ): + return "ElementNull" + + return "Element({0}, {1}, {2})".format( + self.element, self.eltype, self.ghost_type + ) 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/doc/coding_convention/snippet/for b/doc/coding_convention/snippet/for index 8bd190800..4973bc15a 100644 --- a/doc/coding_convention/snippet/for +++ b/doc/coding_convention/snippet/for @@ -1,5 +1,5 @@ #name : for (...; ...; ...) { ... } # -- -for (${1:UInt ${2:i} = 0}; ${3:$2 < N}; ${4:++$2}) { +for (${1:Int ${2:i} = 0}; ${3:$2 < N}; ${4:++$2}) { $0 } diff --git a/doc/dev-doc/manual/getting_started.rst b/doc/dev-doc/manual/getting_started.rst index 90ac8abd7..e62a2927c 100644 --- a/doc/dev-doc/manual/getting_started.rst +++ b/doc/dev-doc/manual/getting_started.rst @@ -1,496 +1,496 @@ Getting Started =============== Building ``Akantu`` -------------------- Dependencies ```````````` In order to compile ``Akantu`` any compiler supporting fully C++14 should work. In addition some libraries are required: - CMake (>= 3.5.1) - Boost (preprocessor and Spirit) - zlib - blas/lapack For the python interface: - Python (>=3 is recommended) - pybind11 (if not present the build system will try to download it) To run parallel simulations: - MPI - Scotch To use the static or implicit dynamic solvers at least one of the following libraries is needed: - MUMPS (since this is usually compiled in static you also need MUMPS dependencies) - PETSc To compile the tests and examples: - Gmsh - google-test (if not present the build system will try to download it) Configuring and compilation ``````````````````````````` ``Akantu`` is a `CMake `_ project, so to configure it, you can either follow the usual way:: > cd akantu > mkdir build > cd build > ccmake .. [ Set the options that you need ] > make > make install Using the python interface -------------------------- You can install ``Akantu`` using pip:: > pip install akantu You can then import the package in a python script as:: import akantu The python API is similar to the C++ one, see :ref:`reference` . If you encouter any problem with the python interface, you are welcome to do a merge request or post an issue on `GitLab `_ . Tutorials with the python interface ``````````````````````````````````` To help getting started, several tutorials using the python interface are available as notebooks with pre-installed version of ``Akantu`` on Binder. The following tutorials are currently available: `Plate whith a hole loaded `_ `Loaded cohesive crack `_ `Making your constitutive law in python `_ Writing a ``main`` function --------------------------- ``Akantu`` first needs to be initialized. The memory management included in the core library handles the correct allocation and de-allocation of vectors, structures and/or objects. Moreover, in parallel computations, the initialization procedure performs the communication setup. This is achieved by the function :cpp:func:`initialize ` that is used as follows:: #include "aka_common.hh" #include "..." using namespace akantu; int main(int argc, char *argv[]) { initialize("input_file.dat", argc, argv); // your code ... } The :cpp:func:`initialize ` function takes the text inpute file and the program parameters which can be parsed by ``Akantu`` in due form (see sect:parser). Obviously it is necessary to include all files needed in main. In this manual all provided code implies the usage of ``akantu`` as namespace. Compiling your simulation ------------------------- The easiest way to compile your simulation is to create a ``cmake`` project by putting all your code in some directory of your choosing. Then, make sure that you have ``cmake`` installed and create a ``CMakeLists.txt`` file. An example of a minimal ``CMakeLists.txt`` file would look like this: .. code-block:: cmake project(my_simu) cmake_minimum_required(VERSION 3.0.0) find_package(Akantu REQUIRED) add_akantu_simulation(my_simu my_simu.cc) Then create a directory called ``build`` and inside it execute ``cmake -DAkantu_DIR= -DCMAKE_BUILD_TYPE=RelWithDebInfo ..``. If you installed ``Akantu`` in a standard directory such as ``/usr/local`` (using ``make install``), you can omit the ``-DAkantu_DIR=`` option. Other why ``path_to_akantu`` is either the folder where you built ``Akantu`` if you did not do a ``make install``, or if you installed ``Akantu`` in ``CMAKE_INSTALL_PREFIX`` it is ``/share/cmake/Akantu``. Once ``cmake`` managed to configure and generate a ``makefile`` you can just do ``make`` .. _loading_mesh: Creating and Loading a Mesh --------------------------- In its current state, ``Akantu`` supports three types of meshes: Gmsh, Abaqus and Diana. Once a :cpp:class:`akantu::Mesh` object is created with a given spatial dimension, it can be filled by reading a mesh input file. The method :cpp:func:`read ` of the class :cpp:class:`Mesh ` infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified. :: - UInt spatial_dimension = 2; + Int spatial_dimension = 2; Mesh mesh(spatial_dimension); // Reading Gmsh files mesh.read("my_gmsh_mesh.msh"); mesh.read("my_gmsh_mesh", _miot_gmsh); The Gmsh reader adds the geometrical and physical tags as mesh data. The -physical values are stored as a :cpp:type:`UInt ` data called +physical values are stored as a :cpp:type:`Int ` data called ``tag_0``, if a string name is provided it is stored as a ``std::string`` data -named ``physical_names``. The geometrical tag is stored as a :cpp:type:`UInt -` data named ``tag_1``. +named ``physical_names``. The geometrical tag is stored as a :cpp:type:`Int +` data named ``tag_1``. Using Arrays ------------ Data in ``Akantu`` can be stored in data containers implemented by the :cpp:class:`akantu::Array` class. In its most basic usage, the :cpp:class:`Array ` class implemented in \akantu is similar to the ``std::vector`` class of the Standard Template Library (STL) for C++. A simple :cpp:class:`Array ` containing a sequence of ``nb_element`` values (of a given type) can be generated with:: Array example_array(nb_element); where ``type`` usually is :cpp:type:`Real `, :cpp:type:`Int `, :cpp:type:`UInt ` or ``bool``. Each value is associated to an index, so that data can be accessed by typing:: auto & val = example_array(index); ``Arrays`` can also contain tuples of values for each index. In that case, the number of components per tuple must be specified at the :cpp:class:`Array ` creation. For example, if we want to create an :cpp:class:`Array ` to store the coordinates (sequences of three values) of ten nodes, the appropriate code is the following:: UInt nb_nodes = 10; UInt spatial_dimension = 3; Array position(nb_nodes, spatial_dimension); In this case the :math:`x` position of the eighth node number will be given by ``position(7, 0)`` (in C++, numbering starts at 0 and not 1). If the number of components for the sequences is not specified, the default value of 1 is used. Here is a list of some basic operations that can be performed on :cpp:class:`Array `: - :cpp:func:`resize(size) ` change the size of the :cpp:class:`Array `. - :cpp:func:`clear ` reset the size of the :cpp:class:`Array ` to zero. (*warning* this changed in > v4.0) - :cpp:func:`set(t) ` set all entries of the :cpp:class:`Array ` to ``t``. - :cpp:func:`copy(const Array & other) ` copy another :cpp:class:`Array ` into the current one. The two :cpp:class:`Arrays ` should have the same number of components. - :cpp:func:`push_back(tuple) ` append a tuple with the correct number of components at the end of the :cpp:class:`Array `. - :cpp:func:`erase(i) ` erase the value at the i-th position. - :cpp:func:`find(value) ` search ``value`` in the current :cpp:class:`Array `. Return position index of the first occurence or -1 if not found. - :cpp:func:`storage() ` Return the address of the allocated memory of the :cpp:class:`Array `. Array iterators --------------- It is very common in ``Akantu`` to loop over arrays to perform a specific treatment. This ranges from geometric calculation on nodal quantities to tensor algebra (in constitutive laws for example). The :cpp:class:`Array ` object has the possibility to request iterators in order to make the writing of loops easier and enhance readability. For instance, a loop over the nodal coordinates can be performed like:: // accessing the nodal coordinates Array // with spatial_dimension components const auto & nodes = mesh.getNodes(); for (const auto & coords : make_view(nodes, spatial_dimension)) { // do what you need .... } In that example, each ``coords`` is a :cpp:class:`Vector\ ` containing geometrical array of size ``spatial_dimension`` and the iteration is conveniently performed by the :cpp:class:`Array ` iterator. The :cpp:class:`Array ` object is intensively used to store second order tensor values. In that case, it should be specified that the returned object type is a matrix when constructing the iterator. This is done when calling the :cpp:func:`make_view `. For instance, assuming that we have a :cpp:class:`Array ` storing stresses, we can loop over the stored tensors by:: for (const auto & stress : make_view(stresses, spatial_dimension, spatial_dimension)) { // stress is of type `const Matrix&` } In that last example, the :cpp:class:`Matrix\ ` objects are ``spatial_dimension`` :math:`\times` ``spatial_dimension`` matrices. The light objects :cpp:class:`Matrix\ ` and :cpp:class:`Vector\ ` can be used and combined to do most common linear algebra. If the number of component is 1, it is possible to use :cpp:func:`make_view ` to this effect. In general, a mesh consists of several kinds of elements. Consequently, the amount of data to be stored can differ for each element type. The straightforward example is the connectivity array, namely the sequences of nodes belonging to each element (linear triangular elements have fewer nodes than, say, rectangular quadratic elements etc.). A particular data structure called :cpp:class:`ElementTypeMapArray\ ` is provided to easily manage this kind of data. It consists of a group of ``Arrays``, each associated to an element type. The following code can retrieve the :cpp:class:`ElementTypeMapArray\ ` which stores the connectivity arrays for a mesh:: const ElementTypeMapArray & connectivities = mesh.getConnectivities(); Then, the specific array associated to a given element type can be obtained by:: const Array & connectivity_triangle = connectivities(_triangle_3); where the first order 3-node triangular element was used in the presented piece of code. Vector & Matrix ``````````````` The :cpp:class:`Array\ ` iterators as presented in the previous section can be shaped as :cpp:class:`Vector\ ` or :cpp:class:`Matrix\ `. This objects represent 1st and 2nd order tensors. As such they come with some functionalities that we will present a bit more into detail in this here. ``Vector`` ''''''''''''' - Accessors: - :cpp:func:`v(i) ` gives the ``i`` -th component of the vector ``v`` - :cpp:func:`v[i] ` gives the ``i`` -th component of the vector ``v`` - :cpp:func:`v.size() ` gives the number of component - Level 1: (results are scalars) - :cpp:func:`v.norm() ` returns the geometrical norm (:math:`L_2`) - :cpp:func:`v.norm\() >` returns the :math:`L_N` norm defined as :math:`\left(\sum_i |v(i)|^N\right)^{1/N}`. N can take any positive integer value. There are also some particular values for the most commonly used norms, ``L_1`` for the Manhattan norm, ``L_2`` for the geometrical norm and ``L_inf`` for the norm infinity. - :cpp:func:`v.dot(x) ` return the dot product of ``v`` and ``x`` - :cpp:func:`v.distance(x) ` return the geometrical norm of :math:`v - x` - Level 2: (results are vectors) - :cpp:func:`v += s `, :cpp:func:`v -= s `, :cpp:func:`v *= s `, :cpp:func:`v /= s ` those are element-wise operators that sum, substract, multiply or divide all the component of ``v`` by the scalar ``s`` - :cpp:func:`v += x `, :cpp:func:`v -= x ` sums or substracts the vector ``x`` to/from ``v`` - :cpp:func:`v.mul(A, x, alpha) ` stores the result of :math:`\alpha \boldsymbol{A} \vec{x}` in ``v``, :math:`\alpha` is equal to 1 by default - :cpp:func:`v.solve(A, b) ` stores the result of the resolution of the system :math:`\boldsymbol{A} \vec{x} = \vec{b}` in ``v`` - :cpp:func:`v.crossProduct(v1, v2) ` computes the cross product of ``v1`` and ``v2`` and stores the result in ``v`` ``Matrix`` ''''''''''''' - Accessors: - :cpp:func:`A(i, j) ` gives the component :math:`A_{ij}` of the matrix ``A`` - :cpp:func:`A(i) ` gives the :math:`i^{th}` column of the matrix as a ``Vector`` - :cpp:func:`A[k] ` gives the :math:`k^{th}` component of the matrix, matrices are stored in a column major way, which means that to access :math:`A_{ij}`, :math:`k = i + j M` - :cpp:func:`A.rows() ` gives the number of rows of ``A`` (:math:`M`) - :cpp:func:`A.cols() ` gives the number of columns of ``A`` (:math:`N`) - :cpp:func:`A.size() ` gives the number of component in the matrix (:math:`M \times N`) - Level 1: (results are scalars) - :cpp:func:`A.norm() ` is equivalent to ``A.norm()`` - :cpp:func:`A.norm\() >` returns the :math:`L_N` norm defined as :math:`\left(\sum_i\sum_j |A(i,j)|^N\right)^{1/N}`. N can take any positive integer value. There are also some particular values for the most commonly used norms, ``L_1`` for the Manhattan norm, ``L_2`` for the geometrical norm and ``L_inf`` for the norm infinity. - :cpp:func:`A.trace() ` return the trace of ``A`` - :cpp:func:`A.det() ` return the determinant of ``A`` - :cpp:func:`A.doubleDot(B) ` return the double dot product of ``A`` and ``B``, :math:`\mat{A}:\mat{B}` - Level 3: (results are matrices) - :cpp:func:`A.eye(s) `, ``Matrix::eye(s)`` fills/creates a matrix with the :math:`s\mat{I}` with :math:`\mat{I}` the identity matrix - :cpp:func:`A.inverse(B) ` stores :math:`\mat{B}^{-1}` in ``A`` - :cpp:func:`A.transpose() ` returns :math:`\mat{A}^{t}` - :cpp:func:`A.outerProduct(v1, v2) ` stores :math:`\vec{v_1} \vec{v_2}^{t}` in ``A`` - :cpp:func:`C.mul\(A, B, alpha) `: stores the result of the product of ``A`` and code{B} time the scalar ``alpha`` in ``C``. ``t_A`` and ``t_B`` are boolean defining if ``A`` and ``B`` should be transposed or not. +----------+----------+--------------+ |``t_A`` |``t_B`` |result | | | | | +----------+----------+--------------+ |false |false |:math:`\mat{C}| | | |= \alpha | | | |\mat{A} | | | |\mat{B}` | | | | | +----------+----------+--------------+ |false |true |:math:`\mat{C}| | | |= \alpha | | | |\mat{A} | | | |\mat{B}^t` | | | | | +----------+----------+--------------+ |true |false |:math:`\mat{C}| | | |= \alpha | | | |\mat{A}^t | | | |\mat{B}` | | | | | +----------+----------+--------------+ |true |true |:math:`\mat{C}| | | |= \alpha | | | |\mat{A}^t | | | |\mat{B}^t` | +----------+----------+--------------+ - :cpp:func:`A.eigs(d, V) ` this method computes the eigenvalues and eigenvectors of ``A`` and store the results in ``d`` and ``V`` such that :math:`d(i) = \lambda_i` and :math:`V(i) = \vec{v_i}` with :math:`\mat{A}\vec{v_i} = \lambda_i\vec{v_i}` and :math:`\lambda_1 > ... > \lambda_i > ... > \lambda_N` .. _sect-common-groups: Mesh ---- Manipulating group of nodes and/or elements ``````````````````````````````````````````` ``Akantu`` provides the possibility to manipulate subgroups of elements and nodes. Any :cpp:class:`ElementGroup ` and/or :cpp:class:`NodeGroup ` must be managed by a :cpp:class:`GroupManager `. Such a manager has the role to associate group objects to names. This is a useful feature, in particular for the application of the boundary conditions, as will be demonstrated in section :ref:`sect-smm-boundary`. To most general group manager is the :cpp:class:`Mesh ` class which inherits from :cpp:class:`GroupManager `. For instance, the following code shows how to request an element group to a mesh: .. code-block:: c++ // request creation of a group of nodes NodeGroup & my_node_group = mesh.createNodeGroup("my_node_group"); // request creation of a group of elements ElementGroup & my_element_group = mesh.createElementGroup("my_element_group"); /* fill and use the groups */ The ``NodeGroup`` object '''''''''''''''''''''''' A group of nodes is stored in :cpp:class:`NodeGroup ` objects. They are quite simple objects which store the indexes of the selected nodes in a :cpp:class:`Array\ `. Nodes are selected by adding them when calling :cpp:func:`add `. For instance you can select nodes having a positive :math:`X` coordinate with the following code: .. code-block:: c++ const auto & nodes = mesh.getNodes(); auto & group = mesh.createNodeGroup("XpositiveNode"); for (auto && data : enumerate(make_view(nodes, spatial_dimension))){ auto node = std::get<0>(data); const auto & position = std::get<1>(data); if (position(0) > 0) group.add(node); } The ``ElementGroup`` object ''''''''''''''''''''''''''' A group of elements is stored in :cpp:class:`ElementGroup ` objects. Since a group can contain elements of various types the :cpp:class:`ElementGroup ` object stores indexes in a :cpp:class:`ElementTypeMapArray\ ` object. Then elements can be added to the group by calling :cpp:func:`add `. For instance, selecting the elements for which the barycenter of the nodes has a positive :math:`X` coordinate can be made with: .. code-block:: c++ auto & group = mesh.createElementGroup("XpositiveElement"); Vector barycenter(spatial_dimension); for_each_element(mesh, [&](auto && element) { mesh.getBarycenter(element, barycenter); if (barycenter(_x) > 0.) { group.add(element); } }); diff --git a/doc/manual/manual-changelog.tex b/doc/manual/manual-changelog.tex index 8e7b937f1..5a752203f 100644 --- a/doc/manual/manual-changelog.tex +++ b/doc/manual/manual-changelog.tex @@ -1,47 +1,47 @@ \chapter*{Changes, New Features, and Fixes} \section*{Version 3.0.0} \begin{itemize} \item[\textbf{\texttt{c++14}}] Switch from C++ standard \code{2003} to \code{2014}\\ Example of changes implied by this: \begin{cpp} - for (UInt g = _not_ghost; g <= _ghost; ++g) { + for (Int g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType)g; Mesh::type_iterator it = this->mesh.firstType(spatial_dimension, gt); Mesh::type_iterator end = this->mesh.lastType(spatial_dimension, gt); for (; it != end; ++it) { ElementType & type = *it; ... } } \end{cpp} becomes \begin{cpp} for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(spatial_dimension, ghost_type)) { ... } } \end{cpp} \item[\textbf{\texttt{feature}}] Parallel cohesive elements \item[\textbf{\texttt{feature}}] Models using new interface for solvers \begin{itemize} \item Same configuration for all models \item Solver can be configured in input file \item PETSc interface temporary inactive \item Periodic boundary condition temporary inactive \end{itemize} \item[\textbf{\texttt{feature}}] Element groups created by default for \code{``physical\_names''} \item[\textbf{\texttt{feature}}] Named arguments for functions (e.g. \code{model.initFull(\_analysis\_method = \_static)}) \item[\textbf{\texttt{api}}] Only one function to solve a step \code{model.solveStep()} \item[\textbf{\texttt{api}}] Simplification of the parallel simulation with the \code{mesh.distribute()} function \end{itemize} %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/doc/manual/manual-feengine.tex b/doc/manual/manual-feengine.tex index ef64c823e..4e5d2649c 100644 --- a/doc/manual/manual-feengine.tex +++ b/doc/manual/manual-feengine.tex @@ -1,75 +1,75 @@ \chapter{FEEngine\index{FEEngine}} \label{chap:feengine} The \code{FEEngine} interface is dedicated to handle the finite-element approximations and the numerical integration of the weak form. As we will see in Chapter \ref{sect:smm}, \code{Model} creates its own \code{FEEngine} object so the explicit creation of the object is not required. \section{Mathematical Operations\label{sect:fe:mathop}} Using the \code{FEEngine} object, one can compute a interpolation, an integration or a gradient. A simple example is given below. \begin{cpp} // having a FEEngine object FEEngine *fem = new FEEngineTemplate(my_mesh, dim, "my_fem"); // instead of this, a FEEngine object can be get using the model: // model.getFEEngine() //compute the gradient Array u; //append the values you want Array nablauq; //gradient array to be computed // compute the gradient fem->gradientOnIntegrationPoints(const Array &u, Array &nablauq, - const UInt nb_degree_of_freedom, + const Int nb_degree_of_freedom, ElementType type); // interpolate Array uq; //interpolated array to be computed // compute the interpolation fem->interpolateOnIntegrationPoints(const Array &u, Array &uq, - UInt nb_degree_of_freedom, + Int nb_degree_of_freedom, ElementType type); // interpolated function can be integrated over the elements Array int_val_on_elem; // integrate fem->integrate(const Array &uq, Array &int_uq, - UInt nb_degree_of_freedom, + Int nb_degree_of_freedom, ElementType type); \end{cpp} Another example below shows how to integrate stress and strain fields over elements assigned to a particular material. \begin{cpp} -UInt sp_dim = 3; //spatial dimension -UInt m = 1; //material index of interest +Int sp_dim = 3; //spatial dimension +Int m = 1; //material index of interest const ElementType type = _tetrahedron_4; //element type // get the stress and strain arrays associated to the material index m const Array & strain_vec = model.getMaterial(m).getGradU(type); const Array & stress_vec = model.getMaterial(m).getStress(type); // get the element filter for the material index const Array & elem_filter = model.getMaterial(m).getElementFilter(type); // initialize the integrated stress and strain arrays Array int_strain_vec(elem_filter.getSize(), sp_dim*sp_dim, "int_of_strain"); Array int_stress_vec(elem_filter.getSize(), sp_dim*sp_dim, "int_of_stress"); // integrate the fields model.getFEEngine().integrate(strain_vec, int_strain_vec, sp_dim*sp_dim, type, _not_ghost, elem_filter); model.getFEEngine().integrate(stress_vec, int_stress_vec, sp_dim*sp_dim, type, _not_ghost, elem_filter); \end{cpp} \input{manual-elements} diff --git a/doc/manual/manual-gettingstarted.tex b/doc/manual/manual-gettingstarted.tex index 573f19639..8d7d873d4 100644 --- a/doc/manual/manual-gettingstarted.tex +++ b/doc/manual/manual-gettingstarted.tex @@ -1,466 +1,466 @@ \chapter{Getting Started} \section{Downloading the Code} The \akantu source code can be requested using the form accessible at the URL \url{http://lsms.epfl.ch/akantu}. There, you will be asked to accept the LGPL license terms. \section{Compiling \akantu} \akantu is a \code{cmake} project, so to configure it, you can either follow the usual way: \begin{command} > cd akantu > mkdir build > cd build > ccmake .. [ Set the options that you need ] > make > make install \end{command} \noindent Or, use the \code{Makefile} we added for your convenience to handle the \code{cmake} configuration \begin{command} > cd akantu > make config > make > make install \end{command} \noindent All the \akantu options are documented in Appendix \ref{app:package-dependencies}. \section{Writing a \texttt{main} Function\label{sect:common:main}} \label{sec:writing_main} First of all, \akantu needs to be initialized. The memory management included in the core library handles the correct allocation and de-allocation of vectors, structures and/or objects. Moreover, in parallel computations, the initialization procedure performs the communication setup. This is achieved by a pair of functions (\code{initialize} and \code{finalize}) that are used as follows: \begin{cpp} #include "aka_common.hh" #include "..." using namespace akantu; int main(int argc, char *argv[]) { initialize("input_file.dat", argc, argv); // your code ... finalize(); } \end{cpp} The \code{initialize} function takes the text inpute file and the program parameters which can be parsed by \akantu in due form (see \ref{sect:parser}). Obviously it is necessary to include all files needed in main. In this manual all provided code implies the usage of \code{akantu} as namespace. \section{Creating and Loading a Mesh\label{sect:common:mesh}} In its current state, \akantu supports three types of meshes: Gmsh~\cite{gmsh}, Abaqus~\cite{abaqus} and Diana~\cite{diana}. Once a \code{Mesh} object is created with a given spatial dimension, it can be filled by reading a mesh input file. The method \code{read} of the class \code{Mesh} infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified. \begin{cpp} -UInt spatial_dimension = 2; +Int spatial_dimension = 2; Mesh mesh(spatial_dimension); // Reading Gmsh files mesh.read("my_gmsh_mesh.msh"); mesh.read("my_gmsh_mesh", _miot_gmsh); // Reading Abaqus files mesh.read("my_abaqus_mesh.inp"); mesh.read("my_abaqus_mesh", _miot_abaqus); // Reading Diana files mesh.read("my_diana_mesh.dat"); mesh.read("my_diana_mesh", _miot_diana); \end{cpp} The Gmsh reader adds the geometrical and physical tags as mesh data. The physical -values are stored as a \code{UInt} data called \code{tag\_0}, if a string +values are stored as a \code{Int} data called \code{tag\_0}, if a string name is provided it is stored as a \code{std::string} data named -\code{physical\_names}. The geometrical tag is stored as a \code{UInt} data named +\code{physical\_names}. The geometrical tag is stored as a \code{Int} data named \code{tag\_1}. The Abaqus reader stores the \code{ELSET} in ElementGroups and the \code{NSET} in NodeGroups. The material assignment can be retrieved from the \code{std::string} mesh data named \code{abaqus\_material}. % \akantu supports meshes generated with Gmsh~\cite{gmsh}, a free % software available at \url{http://geuz.org/gmsh/} where a detailed % documentation can be found. Consequently, this manual will not provide % Gmsh usage directions. Gmsh outputs meshes in \code{.msh} format that can be read % by \akantu. In order to import a mesh, it is necessary to create % a \code{Mesh} object through the following function calls: % \begin{cpp} -% UInt spatial_dimension = 2; +% Int spatial_dimension = 2; % Mesh mesh(spatial_dimension); % \end{cpp} % The only parameter that has to be specified by the user is the spatial % dimension of the problem. Now it is possible to read a \code{.msh} file with % a \code{MeshIOMSH} object that takes care of loading a mesh to memory. % This step is carried out by: % \begin{cpp} % mesh.read("square.msh"); % \end{cpp} % where the \code{MeshIOMSH} object is first created before being % used to read the \code{.msh} file. The mesh file name as well as the \code{Mesh} % object must be specified by the user. % The \code{MeshIOMSH} object can also write mesh files. This feature % is useful to save a mesh that has been modified during a % simulation. The \code{write} method takes care of it: % \begin{cpp} % mesh_io.write("square_modified.msh", mesh); % \end{cpp} % which works exactly like the \code{read} method. % \akantu supports also meshes generated by % DIANA (\url{http://tnodiana.com}), but only in reading mode. A similar % procedure applies where the only % difference is that the \code{MeshIODiana} object should be used % instead of the \code{MeshIOMSH} one. Additional mesh readers can be % introduced into \akantu by coding new \code{MeshIO} classes. \section{Using \texttt{Arrays}} Data in \akantu can be stored in data containers implemented by the \code{Array} class. In its most basic usage, the \code{Array} class implemented in \akantu is similar to the \code{vector} class of the Standard Template Library (STL) for C++. A simple \code{Array} containing a sequence of \code{nb\_element} values (of a given type) can be generated with: \begin{cpp} Array example_array(nb_element); \end{cpp} where \code{type} usually is \code{Real}, \code{Int}, \code{UInt} or \code{bool}. Each value is associated to an index, so that data can be accessed by typing: \begin{cpp} auto & val = example_array(index) \end{cpp} \code{Arrays} can also contain tuples of values for each index. In that case, the number of components per tuple must be specified at the \code{Array} creation. For example, if we want to create an \code{Array} to store the coordinates (sequences of three values) of ten nodes, the appropriate code is the following: \begin{cpp} - UInt nb_nodes = 10; - UInt spatial_dimension = 3; + Int nb_nodes = 10; + Int spatial_dimension = 3; Array position(nb_nodes, spatial_dimension); \end{cpp} In this case the $x$ position of the eighth node number will be given by \code{position(7, 0)} (in C++, numbering starts at 0 and not 1). If the number of components for the sequences is not specified, the default value of 1 is used. Here is a list of some basic operations that can be performed on \code{Array}: \begin{itemize} \item \code{resize(size)} change the size of the \code{Array}. \item \code{clear()} set all entries of the \code{Array} to zero. \item \code{set(t)} set all entries of the \code{Array} to \code{t}. \item \code{copy(const Array $\&$ other)} copy another \code{Array} into the current one. The two \code{Array} should have the same number of components. \item \code{push$\_$back(tuple)} append a tuple with the correct number of components at the end of the \code{Array}. \item \code{erase(i) erase the value at the i-th position.} \item \code{find(value)} search \code{value} in the current \code{Array}. Return position index of the first occurence or $-1$ if not found. \item \code{storage()} Return the address of the allocated memory of the \code{Array}. \end{itemize} \subsection{\texttt{Arrays} iterators} It is very common in \akantu to loop over arrays to perform a specific treatment. This ranges from geometric calculation on nodal quantities to tensor algebra (in constitutive laws for example). The \code{Array} object has the possibility to request iterators in order to make the writing of loops easier and enhance readability. For instance, a loop over the nodal coordinates can be performed like: \begin{cpp} //accessing the nodal coordinates Array (spatial_dimension components) const auto & nodes = mesh.getNodes(); //creating the iterators auto it = nodes.begin(spatial_dimension); auto end = nodes.end(spatial_dimension); for (; it != end; ++it){ const auto & coords = (*it); //do what you need .... } \end{cpp} In that example, each \code{coords} is a \code{Vector} containing geometrical array of size \code{spatial\_dimension} and the iteration is conveniently performed by the \code{Array} iterator. With the switch to \code{c++14} this can be also written as \begin{cpp} //accessing the nodal coordinates Array (spatial_dimension components) const auto & nodes = mesh.getNodes(); for (const auto & coords : make_view(nodes, spatial_dimension) { //do what you need .... } \end{cpp} The \code{Array} object is intensively used to store second order tensor values. In that case, it should be specified that the returned object type is a matrix when constructing the iterator. This is done when calling the \code{begin} function. For instance, assuming that we have a \code{Array} storing stresses, we can loop over the stored tensors by: \begin{cpp} //creating the iterators auto it = stresses.begin(spatial_dimension, spatial_dimension); auto end = stresses.end(spatial_dimension, spatial_dimension); for (; it != end; ++it){ Matrix & stress = (*it); //do what you need .... } \end{cpp} In that last example, the \code{Matrix} objects are \code{spatial\_dimension} $\times$ \code{spatial\_dimension} matrices. The light objects \code{Matrix} and \code{Vector} can be used and combined to do most common linear algebra. If the number of component is 1, it is possible to use a scalar\_iterator rather than the vector/matrix one. In general, a mesh consists of several kinds of elements. Consequently, the amount of data to be stored can differ for each element type. The straightforward example is the connectivity array, namely the sequences of nodes belonging to each element (linear triangular elements have fewer nodes than, say, rectangular quadratic elements etc.). A particular data structure called \code{ElementTypeMapArray} is provided to easily manage this kind of data. It consists of a group of \code{Arrays}, each associated to an element type. The following code can retrieve the \code{ElementTypeMapArray} which stores the connectivity arrays for a mesh: \begin{cpp} const ElementTypeMapArray & connectivities = mesh.getConnectivities(); \end{cpp} Then, the specific array associated to a given element type can be obtained by \begin{cpp} const Array & connectivity_triangle = connectivities(_triangle_3); \end{cpp} where the first order 3-node triangular element was used in the presented piece of code. \subsection{Vector \& Matrix} The \code{Array} iterators as presented in the previous section can be shaped as \code{Vector} or \code{Matrix}. This objects represent $1^{st}$ and $2^{nd}$ order tensors. As such they come with some functionalities that we will present a bit more into detail in this here. \subsubsection{\texttt{Vector}} \begin{enumerate} \item Accessors: \begin{itemize} \item \code{v(i)} gives the $i^{th}$ component of the vector \code{v} \item \code{v[i]} gives the $i^{th}$ component of the vector \code{v} \item \code{v.size()} gives the number of component \end{itemize} \item Level 1: (results are scalars) \begin{itemize} \item \code{v.norm()} returns the geometrical norm ($L_2$) \item \code{v.norm()} returns the $L_N$ norm defined as $\left(\sum_i |\code{v(i)}|^N\right)^{1/N}$. N can take any positive integer value. There are also some particular values for the most commonly used norms, \code{L\_1} for the Manhattan norm, \code{L\_2} for the geometrical norm and \code{L\_inf} for the norm infinity. \item \code{v.dot(x)} return the dot product of \code{v} and \code{x} \item \code{v.distance(x)} return the geometrical norm of $\code{v} - \code{x}$ \end{itemize} \item Level 2: (results are vectors) \begin{itemize} \item \code{v += s}, \code{v -= s}, \code{v *= s}, \code{v /= s} those are element-wise operators that sum, substract, multiply or divide all the component of \code{v} by the scalar \code{s} \item \code{v += x}, \code{v -= x} sums or substracts the vector \code{x} to/from \code{v} \item \code{v.mul(A, x, alpha)} stores the result of $\alpha \mat{A} \vec{x}$ in \code{v}, $\alpha$ is equal to 1 by default \item \code{v.solve(A, b)} stores the result of the resolution of the system $\mat{A} \vec{x} = \vec{b}$ in \code{v} \item \code{v.crossProduct(v1, v2)} computes the cross product of \code{v1} and \code{v2} and stores the result in \code{v} \end{itemize} \end{enumerate} \subsubsection{\texttt{Matrix}} \begin{enumerate} \item Accessors: \begin{itemize} \item \code{A(i, j)} gives the component $A_{ij}$ of the matrix \code{A} \item \code{A(i)} gives the $i^{th}$ column of the matrix as a \code{Vector} \item \code{A[k]} gives the $k^{th}$ component of the matrix, matrices are stored in a column major way, which means that to access $A_{ij}$, $k = i + j M$ \item \code{A.rows()} gives the number of rows of \code{A} ($M$) \item \code{A.cols()} gives the number of columns of \code{A} ($N$) \item \code{A.size()} gives the number of component in the matrix ($M \times N$) \end{itemize} \item Level 1: (results are scalars) \begin{itemize} \item \code{A.norm()} is equivalent to \code{A.norm()} \item \code{A.norm()} returns the $L_N$ norm defined as $\left(\sum_i\sum_j |\code{A(i,j)}|^N\right)^{1/N}$. N can take any positive integer value. There are also some particular values for the most commonly used norms, \code{L\_1} for the Manhattan norm, \code{L\_2} for the geometrical norm and \code{L\_inf} for the norm infinity. \item \code{A.trace()} return the trace of \code{A} \item \code{A.det()} return the determinant of \code{A} \item \code{A.doubleDot(B)} return the double dot product of \code{A} and \code{B}, $\mat{A}:\mat{B}$ \end{itemize} \item Level 3: (results are matrices) \begin{itemize} \item \code{A.eye(s)}, \code{Matrix::eye(s)} fills/creates a matrix with the $s\mat{I}$ with $\mat{I}$ the identity matrix \item \code{A.inverse(B)} stores $\mat{B}^{-1}$ in \code{A} \item \code{A.transpose()} returns $\mat{A}^{t}$ \item \code{A.outerProduct(v1, v2)} stores $\vec{v_1} \vec{v_2}^{t}$ in \code{A} \item \code{C.mul(A, B, alpha)}: stores the result of the product of \code{A} and code{B} time the scalar \code{alpha} in \code{C}. \code{t\_A} and \code{t\_B} are boolean defining if \code{A} and \code{B} should be transposed or not. \begin{tabular}{ccl} \toprule \code{t\_A} & \code{t\_B} & result \\ \midrule false & false & $\mat{C} = \alpha \mat{A} \mat{B}$\\ false & true & $\mat{C} = \alpha \mat{A} \mat{B}^t$\\ true & false & $\mat{C} = \alpha \mat{A}^t \mat{B}$\\ true & true & $\mat{C} = \alpha \mat{A}^t \mat{B}^t$\\ \bottomrule \end{tabular} \item \code{A.eigs(d, V)} this method computes the eigenvalues and eigenvectors of \code{A} and store the results in \code{d} and \code{V} such that $\code{d(i)} = \lambda_i$ and $\code{V(i)} = \vec{v_i}$ with $\mat{A}\vec{v_i} = \lambda_i\vec{v_i}$ and $\lambda_1 > ... > \lambda_i > ... > \lambda_N$ \end{itemize} \end{enumerate} \subsubsection{\texttt{Tensor3}} Accessors: \begin{itemize} \item \code{t(i, j, k)} gives the component $T_{ijk}$ of the tensor \code{t} \item \code{t(k)} gives the $k^{th}$ two-dimensional tensor as a \code{Matrix} \item \code{t[k]} gives the $k^{th}$ two-dimensional tensor as a \code{Matrix} \end{itemize} \section{Manipulating group of nodes and/or elements\label{sect:common:groups}} \akantu provides the possibility to manipulate subgroups of elements and nodes. Any \code{ElementGroup} and/or \code{NodeGroup} must be managed by a \code{GroupManager}. Such a manager has the role to associate group objects to names. This is a useful feature, in particular for the application of the boundary conditions, as will be demonstrated in section \ref{sect:smm:boundary}. To most general group manager is the \code{Mesh} class which inheritates from the \code{GroupManager} class. For instance, the following code shows how to request an element group to a mesh: \begin{cpp} // request creation of a group of nodes NodeGroup & my_node_group = mesh.createNodeGroup("my_node_group"); // request creation of a group of elements ElementGroup & my_element_group = mesh.createElementGroup("my_element_group"); /* fill and use the groups */ \end{cpp} \subsection{The \texttt{NodeGroup} object} A group of nodes is stored in \code{NodeGroup} objects. They are quite simple objects which store the indexes of the selected nodes in a \code{Array}. Nodes are selected by adding them when calling \code{NodeGroup::add}. For instance you can select nodes having a positive $X$ coordinate with the following code: \begin{cpp} const auto & nodes = mesh.getNodes(); auto & group = mesh.createNodeGroup("XpositiveNode"); for (auto && data : enumerate(make_view(nodes, spatial_dimension))){ auto node = std::get<0>(data); const auto & position = std::get<1>(data); if (position(0) > 0) group.add(node); } \end{cpp} \subsection{The \texttt{ElementGroup} object} A group of elements is stored in \code{ElementGroup} objects. Since a group can contain elements of various types the \code{ElementGroup} object stores indexes in a \code{ElementTypeMapArray} object. Then elements can be added to the group by calling \code{addElement}. For instance, selecting the elements for which the barycenter of the nodes has a positive $X$ coordinate can be made with: \begin{cpp} auto & group = mesh.createElementGroup("XpositiveElement"); Vector barycenter(spatial_dimension); for(auto type : mesh.elementTypes()){ - UInt nb_element = mesh.getNbElement(type); + Int nb_element = mesh.getNbElement(type); - for(UInt e = 0; e < nb_element; ++e) { + for(Int e = 0; e < nb_element; ++e) { Element element{type, e, _not_ghost}; mesh.getBarycenter(element, barycenter); if (barycenter(_x) > 0.) group.add(element); } } \end{cpp} \section{Compiling your simulation} The easiest way to compile your simulation is to create a \code{cmake} project by putting all your code in some directory of your choosing. Then, make sure that you have \code{cmake} installed and create a \code{CMakeLists.txt} file. An example of a minimal \code{CMakeLists.txt} file would look like this: \begin{cmake} project(my_simu) cmake_minimum_required(VERSION 3.0.0) find_package(Akantu REQUIRED) add_akantu_simulation(my_simu my_simu.cc) \end{cmake} % Then create a directory called \code{build} and inside it execute \code{ccmake ..} which opens a configuration screen. If you installed \akantu in a standard directory such as \code{/usr/local} (using \code{make install}), you should be able to compile by hitting the \code{c} key, setting \code{CMAKE\textunderscore{}BUILD\textunderscore{}TYPE} to \code{Release}, hitting the \code{c} key again a few times and then finishing the configuration with the \code{g} key. After that, you can type \code{make} to build your simulation. If you change your simulation code later on, you only need to type \code{make} again. If you get an error that \code{FindAkantu.cmake} was not found, you have to set the \code{Akantu\textunderscore{}DIR} variable, which will appear after dismissing the error message. If you built \akantu without running \code{make install}, the variable should be set to the \code{build} subdirectory of the \akantu source code. If you installed it in \code{\$PREFIX}, set the variable to \code{\$PREFIX/share/cmake/Akantu} instead. %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/doc/manual/manual-heattransfermodel.tex b/doc/manual/manual-heattransfermodel.tex index cd980f706..d95a6d339 100644 --- a/doc/manual/manual-heattransfermodel.tex +++ b/doc/manual/manual-heattransfermodel.tex @@ -1,163 +1,163 @@ \chapter{Heat Transfer Model} The heat transfer model is a specific implementation of the \code{Model} interface dedicated to handle the dynamic heat equation. \section{Theory} The strong form of the dynamic heat equation can be expressed as \begin{equation} \rho c_v \dot{T} + \nabla \cdot \vec{\kappa} \nabla T = b \end{equation} with $T$ the scalar temperature field, $c_v$ the specific heat capacity, $\rho$ the mass density, $\mat{\kappa}$ the conductivity tensor, and $b$ the heat generation per unit of volume. The discretized weak form with a finite number of elements is \begin{equation} \forall i \quad \sum_j \left( \int_\Omega \rho c_v N_j N_i d\Omega \right) \dot{T}_j - \sum_j \left( \int_\Omega \vec{\kappa} \nabla N_j \nabla N_i d\Omega \right) T_j = - \int_{\Gamma} N_i \vec{q} \cdot \vec{n} d\Gamma + \int_\Omega b N_i d\Omega \end{equation} with $i$ and $j$ the node indices, $\vec{n}$ the normal field to the surface $\Gamma = \partial \Omega$. To simplify, we can define the capacity and the conductivity matrices as \begin{equation} C_{ij} = \int_\Omega \rho c_v N_j N_i d\Omega \qquad \textrm{and} \qquad K_{ij} = - \int_\Omega \vec{\kappa} \nabla N_j \nabla N_i d\Omega \end{equation} and the system to solve can be written \begin{equation} \mat{C} \cdot \vec{\dot{T}} = \vec{Q}^{\text{ext}} -\mat{K} \cdot \vec{T}~, \end{equation} with $\vec{Q}^{\text{ext}}$ the consistent heat generated. \section{Using the Heat Transfer Model} A material file name has to be provided during initialization. Currently, the \code{HeatTransferModel} object uses dynamic analysis with an explicit time integration scheme. It can simply be created like this \begin{cpp} HeatTransferModel model(mesh, spatial_dimension); \end{cpp} while an existing mesh has been used (see \ref{sect:common:mesh}). Then the model object can be initialized with: \begin{cpp} model.initFull() \end{cpp} This function will load the material properties, and allocate / initialize the nodes and element \code{Array}s More precisely, the heat transfer model contains 4 \code{Arrays}: \begin{description} \item[temperature] contains the nodal temperature $T$ (zero by default after the initialization). \item[temperature\_rate] contains the variations of temperature $\dot{T}$ (zero by default after the initialization). \item[blocked\_dofs] contains a Boolean value for each degree of freedom specifying whether the degree is blocked or not. A Dirichlet boundary condition ($T_d$) can be prescribed by setting the \textbf{blocked\_dofs} value of a degree of freedom to \code{true}. The \textbf{temperature} and the \textbf{temperature\_rate} are computed for all degrees of freedom where the \textbf{blocked\_dofs} value is set to \code{false}. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept. \item[external\_heat\_rate] contains the external heat generations. $\vec{Q^{ext}}$ on the nodes. \item[internal\_heat\_rate] contains the internal heat generations. $\vec{Q^{int}} = -\mat{K} \cdot \vec{T}$ on the nodes. \end{description} Only a single material can be specified on the domain. A material text file (\eg material.dat) provides the material properties as follows: \begin{cpp} model heat_transfer_model [ capacity = %\emph{XXX}% density = %\emph{XXX}% conductivity = [%\emph{XXX}% ... %\emph{XXX}%] ] \end{cpp} where the \code{capacity} and \code{density} are scalars, and the \code{conductivity} is specified as a $3\times 3$ tensor. \subsection{Explicit Dynamic} The explicit time integration scheme in \akantu uses a lumped capacity matrix $\mat{C}$ (reducing the computational cost, see Chapter \ref{sect:smm}). This matrix is assembled by distributing the capacity of each element onto its nodes. Therefore, the resulting $\mat{C}$ is a diagonal matrix stored in the \code{capacity} \code{Array} of the model. \begin{cpp} model.assembleCapacityLumped(); \end{cpp} \index{HeatTransferModel!assembleCapacityLumped} \note{Currently, only the explicit time integration with lumped capacity matrix is implemented within \akantu.} The explicit integration scheme is \index{Forward Euler} \emph{Forward Euler} \cite{curnier92a}. \begin{itemize} \item Predictor: $\vec{T}_{n+1} = \vec{T}_{n} + \Delta t \dot{\vec{T}}_{n}$ \item Update residual: $\vec{R}_{n+1} = \left( \vec{Q^{ext}_{n+1}} - \vec{K}\vec{T}_{n+1} \right)$ \item Corrector : $\dot{\vec{T}}_{n+1} = \mat{C}^{-1} \vec{R}_{n+1}$ \end{itemize} The explicit integration scheme is conditionally stable. The time step has to be smaller than the stable time step, and it can be obtained in \akantu as follows: \begin{cpp} time_step = model.getStableTimeStep(); \end{cpp} \index{HeatTransferModel!StableTimeStep} The stable time step is defined as: \begin{equation}\label{eqn:htm:explicit:stabletime} \Delta t_{\st{crit}} = 2 \Delta x^2 \frac{\rho c_v}{\mid\mid \mat{\kappa} \mid\mid^\infty} \end{equation} where $\Delta x$ is the characteristic length (\eg the inradius in the case of linear triangle element), $\rho$ is the density, $\mat{\kappa}$ is the conductivity tensor, and $c_v$ is the specific heat capacity. It is necessary to impose a time step which is smaller than the stable time step, for instance, by multiplying the stable time step by a safety factor smaller than one. \begin{cpp} const Real safety_time_factor = 0.1; Real applied_time_step = time_step * safety_time_factor; model.setTimeStep(applied_time_step); \end{cpp} The following loop allows, for each time step, to update the \code{temperature}, \code{residual} and \code{temperature\_rate} fields following the previously described integration scheme. \begin{cpp} - for (UInt s = 1; (s-1)*applied_time_step < total_time; ++s) { + for (Int s = 1; (s-1)*applied_time_step < total_time; ++s) { model.solveStep(); } \end{cpp} \index{HeatTransferModel!solveStep()} An example of explicit dynamic heat propagation is presented in \\ \shellcode{\examplesdir/heat\_transfer/explicit\_heat\_transfer.cc}. \\ This example consists of a square 2D plate of \SI{1}{\metre^2} having an initial temperature of \SI{100}{\kelvin} everywhere but a none centered hot point maintained at \SI{300}{\kelvin}. Figure~\ref{fig:htm:explicit:dynamic} presents the geometry of this case. The material used is a linear fictitious elastic material with a density of \SI{8940}{\kilo\gram\per\metre^3}, a conductivity of \SI{401}{\watt\per\metre\per\kelvin} and a specific heat capacity of \SI{385}{\joule\per\kelvin\per\kilogram}. The time step used is \SI{0.12}{\second}. \begin{figure}[!htb] \centering \includegraphics[width=.4\textwidth]{figures/hot-point-1} \hfill \includegraphics[width=.4\textwidth]{figures/hot-point-2} \caption{Initial temperature field (left) and after 15000 time steps = 30 minutes (right). The lines represent iso-surfaces.} \label{fig:htm:explicit:dynamic} \end{figure} %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/doc/manual/manual-io.tex b/doc/manual/manual-io.tex index b4aac0022..f0d5284a5 100644 --- a/doc/manual/manual-io.tex +++ b/doc/manual/manual-io.tex @@ -1,297 +1,297 @@ \chapter{Input/Output}\index{I\/O} \section{Input file \label{sect:parser}} The text input file of a simulation should be precised using the method \code{initialize} which will instantiate the static \code{Parser} object of \akantu. This section explains how to manipulate \code{Parser} objects to input data in \akantu. \begin{cpp} int main(int argc, char *argv[]) { initialize("input_files.dat", argc, argv); ... \end{cpp} \subsection{Akantu Parser} \akantu file parser has a tree organization. \begin{itemize} \item \code{Parser}, the root of the tree, can be accessed using \begin{cpp} Parser & parser = getStaticParser(); \end{cpp} \item \code{ParserSection}, branch of the tree, contains map a of sub-sections (\code{SectionType}, \code{ParserSection}) and a \code{ParserSection *} pointing to the parent section. The user section of the input file can directly be accessed by \begin{cpp} const ParserSection & usersect = getUserParser(); \end{cpp} \item \code{ParserParameter}, the leaf of the tree, carries data of the input file which can be casted to the correct type with \begin{cpp} Real mass = usersect.getParameter("mass"); \end{cpp} or used directly within an expression \begin{cpp} Real dead_weight = 9.81 * usersect.getParameterValue("mass"); \end{cpp} \end{itemize} \subsection{Grammar} The structure of text input files consists of different sections containing a list of parameters. As example, the file parsed in the previous section will look like \begin{cpp} user parameters [ mass = 10.5 ] \end{cpp} Basically every standard arithmetic operations can be used inside of input files as well as the constant \code{pi} and \code{e} and the exponent operator \code{\^{}}. Operations between \code{ParserParameter} are also possible with the convention that only parameters of the current and the parent sections are available. \code{Vector} and \code{Matrix} can also be read according to the \code{NumPy}\cite{numpy} writing convention (a.e. cauchy$\_$stress$\_$tensor = [[$\sigma_{xx}$, $\sigma_{xy}$],[$\sigma_{yx}$,$\sigma_{yy}$]]). An example illustrating how to parse the following input file can be found in \code{example$\backslash$io$\backslash$parser$\backslash$example$\_$parser.cc}. \begin{cpp} user parameters [ spatial$\_$dimension = 2 mesh$\_$file = swiss$\_$cheese.msh inner$\_$holes = holes outter$\_$crust = crust lactostatic$\_$p = 30e3 stress = [[lactostatic$\_$p,0],[0,lactostatic$\_$p]] max$\_$nb$\_$iterations = 100 precision = 1e-9 ] \end{cpp} \subsection{Material section \label{sect:io:material}} The input file should also be used to specify material characteristics (constitutive behavior and material properties). The dedicated material section is then read by \code{initFull} method of \code{SolidMechanicsModel} which initializes the different materials specified with the following convention: \begin{cpp} material %\emph{constitutive\_law}% %\emph{}% [ name = $value$ rho = $value$ ... ] \end{cpp} \index{Constitutive\_laws} where \emph{constitutive\_law} is the adopted constitutive law, followed by the material properties listed one by line in the bracket (\eg \code{name} and density \code{rho}). Some constitutive laws can also have an \emph{optional flavor}. More information can be found in sections relative to material constitutive laws \ref{sect:smm:CL} or in Appendix \ref{app:material-parameters}. \section{Output data} \subsection{Generic data} In this chapter, we address ways to get the internal data in human-readable formats. The models in \akantu handle data associated to the mesh, but this data can be split into several \code{Arrays}. For example, the data stored per element type in a \code{ElementTypeMapArray} is composed of as many \code{Array}s as types in the mesh. In order to get this data in a visualization software, the models contain a object to dump \code{VTK} files. These files can be visualized in software such as \code{ParaView}\cite{paraview}, \code{ViSit}\cite{visit} or \code{Mayavi}\cite{mayavi}. The internal dumper of the model can be configured to specify which data fields are to be written. This is done with the \code{addDumpField}\index{I\/O!addDumpField} method. By default all the files are generated in a folder called \code{paraview/} \begin{cpp} model.setBaseName("output"); // prefix for all generated files model.addDumpField("displacement"); model.addDumpField("stress"); ... model.dump() \end{cpp} The fields are dumped with the number of components of the memory. For example, in 2D, the memory has \code{Vector}s of 2 components, or the $2^{nd}$ order tensors with $2\times2$ components. This memory can be dealt with \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} which always dumps \code{Vector}s with 3 components or \code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} which dumps $2^{nd}$ order tensors with $3\times3$ components respectively. The routines \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} and \code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} were introduced because of \code{ParaView} which mostly manipulate 3D data. Those fields which are stored by quadrature point are modified to be seen in the \code{VTK} file as elemental data. To do this, the default is to average the values of all the quadrature points. The list of fields depends on the models (for \code{SolidMechanicsModel} see table~\ref{tab:io:smm_field_list}). \begin{table} \centering \begin{tabular}{llll} \toprule key & type & support \\ \midrule displacement & Vector & nodes \\ mass & Vector & nodes \\ velocity & Vector & nodes \\ acceleration & Vector & nodes \\ force & Vector & nodes \\ residual & Vector & nodes \\ increment & Vector & nodes \\ {blocked\_dofs} & Vector & nodes \\ partitions & Real & elements \\ material\_index & variable & elements\\ strain & Matrix & quadrature points \\ Green strain & Matrix & quadrature points \\ principal strain & Vector & quadrature points \\ principal Green strain & Vector & quadrature points \\ grad\_u & Matrix & quadrature points \\ stress & Matrix & quadrature points \\ Von Mises stress & Real & quadrature points \\ material\_index & variable & quadrature points \\ \bottomrule \end{tabular} \caption{List of dumpable fields for \code{SolidMechanicsModel}.} \label{tab:io:smm_field_list} \end{table} \subsection{Cohesive elements' data} Cohesive elements and their relative data can be easily dumped thanks to a specific dumper contained in \code{SolidMechanicsModelCohesive}. In order to use it, one has just to add the string \code{"cohesive elements"} when calling each method already illustrated. Here is an example on how to dump displacement and damage: \begin{cpp} model.setBaseNameToDumper("cohesive elements", "cohesive_elements_output"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); ... model.dump("cohesive elements"); \end{cpp} \subsubsection{Fragmentation data} Whenever the \code{SolidMechanicsModelCohesive} is used, it is possible to dump additional data about the fragments that get formed in the simulation both in serial and parallel. This task is carried out by the \code{FragmentManager} class, that takes care of computing the following quantities for each fragment: \begin{itemize} \item index; \item mass; \item moments of inertia; \item velocity; \item number of elements. \end{itemize} These computations can be realized at once by calling the function \code{computeAllData}, or individually by calling the other public functions of the class. The data can be dumped to be visualized in Paraview, or can be accessed within the simulation. An example of usage is: \begin{cpp} FragmentManager fragment_manager(model); fragment_manager.buildAllData(); ... model.addDumpField("fragments"); // this field contains the indices model.addDumpField("fragments mass"); model.addDumpField("moments of inertia"); model.addDumpField("elements per fragment"); ... - for (UInt step = 1; step <= total_steps; ++step) { + for (Int step = 1; step <= total_steps; ++step) { ... fragment_manager.buildAllData(); model.dump(); } ... const Array & fragment_velocities = fragment_manager.getVelocity(); ... \end{cpp} At the end of this example the velocities of the fragments are accessed with a reference to a \code{const Array}. The size of this array is the number of fragments, and its number of components is the spatial dimension in this case. \subsection{Advanced dumping} \subsubsection{Arbitrary fields} In addition to the predetermined fields from the models and materials, the user can add any data to a dumper as long as the support is the same. That is to say data that have the size of the full mesh on if the dumper is dumping the mesh, or of the size of an element group if it is a filtered dumper. For this the easiest is to use the ``external'' fields register functions\index{I\/O!addDumpFieldExternal} The simple case force nodal and elemental data are to pass directly the data container itself if it as the good size. \begin{itemize} \item For nodal fields : \begin{cpp} template addDumpFieldExternal(const std::string & field_name, const Array & field); \end{cpp} It is assumed that the array as the same size as the number of nodes in the mesh \begin{cpp} Array nodal_data(nb_nodes, nb_component); mesh.addDumpFieldExternal("my_field", nodal_data); \end{cpp} \item For elemental fields : \begin{cpp} template addDumpFieldExternal(const std::string & field_name, const ElementTypeMapArray & field); \end{cpp} It is assumed that the arrays in the map have the same sizes as the element numbers in the mesh for element types of dimension \code{spatial\_dimension}. \begin{cpp} ElementTypeMapArray elem_data; elem_data.initialize(mesh, _nb_component = 1, _with_nb_element = true); mesh.addDumpFieldExternal("my_field", elem_data); \end{cpp} \end{itemize} If some changes have to be applied on the data as for example a padding for \code{ParaView} vectors, this can be done by using the field interface. \begin{cpp} mesh.addDumpFieldExternal(const std::string & field_name, std::shared & field); \end{cpp} All these functions use the default dumper registered in the mesh but also have the \code{ToDumper} variation with the dumper name specified. For example: \begin{cpp} mesh.addDumpFieldExternalToDumper("cohesive elements", "my_field", my_field); \end{cpp} An example of code presenting this interface is present in the \shellcode{\examplesdir/io/dumper}. This interface is part of the \code{Dumpable} class from which the \code{Mesh} inherits. \subsubsection{Creating a new dumper} You can also create you own dumpers, \akantu uses a third-party library in order to write the output files, \code{IOHelper}. \akantu supports the \code{ParaView} format and a Text format defined by \code{IOHelper}. This two files format are handled by the classes \code{DumperParaview}\index{I\/O!DumperParaview} and \code{DumperText}\index{I\/O!DumperText}. In order to use them you can instantiate on of this object in your code. This dumper have a simple interface. You can register a mesh \code{registerMesh}\index{I\/O!registerMesh}, \code{registerFilteredMesh}\index{I\/O!registerFilteredMesh} or a field, \code{registerField}\index{I\/O!registerField}. An example of code presenting this low level interface is present in the \shellcode{\examplesdir/io/dumper}. The different types of \code{Field} that can be created are present in the source folder \shellcode{src/io/dumper}. %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/doc/manual/manual-solidmechanicsmodel.tex b/doc/manual/manual-solidmechanicsmodel.tex index f0fe8760f..5e321d550 100644 --- a/doc/manual/manual-solidmechanicsmodel.tex +++ b/doc/manual/manual-solidmechanicsmodel.tex @@ -1,1138 +1,1138 @@ \chapter{Solid Mechanics Model\index{SolidMechanicsModel}\label{sect:smm}} The solid mechanics model is a specific implementation of the \code{Model} interface dedicated to handle the equations of motion or equations of equilibrium. The model is created for a given mesh. It will create its own \code{FEEngine} object to compute the interpolation, gradient, integration and assembly operations. A \code{SolidMechanicsModel} object can simply be created like this: \begin{cpp} SolidMechanicsModel model(mesh); \end{cpp} where \code{mesh} is the mesh for which the equations are to be solved. A second parameter called \code{spatial\_dimension} can be added after \code{mesh} if the spatial dimension of the problem is different than that of the mesh. This model contains at least the following six \code{Arrays}: \begin{description} \item[blocked\_dofs] contains a Boolean value for each degree of freedom specifying whether that degree is blocked or not. A Dirichlet boundary condition can be prescribed by setting the \textbf{blocked\_dofs} value of a degree of freedom to \code{true}. A Neumann boundary condition can be applied by setting the \textbf{blocked\_dofs} value of a degree of freedom to \code{false}. The \textbf{displacement}, \textbf{velocity} and \textbf{acceleration} are computed for all degrees of freedom for which the \textbf{blocked\_dofs} value is set to \code{false}. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept. \item[displacement] contains the displacements of all degrees of freedom. It can be either a computed displacement for free degrees of freedom or an imposed displacement in case of blocked ones ($\vec{u}$ in the following). \item[velocity] contains the velocities of all degrees of freedom. As \textbf{displacement}, it contains computed or imposed velocities depending on the nature of the degrees of freedom ($\dot{\vec{u}}$ in the following). \item[acceleration] contains the accelerations of all degrees of freedom. As \textbf{displacement}, it contains computed or imposed accelerations depending on the nature of the degrees of freedom ($\ddot{\vec{u}}$ in the following). \item[external\_force] contains the external forces applied on the nodes ($\vec{f}_{\st{ext}}$ in the following). \item[internal\_force] contains the internal forces on the nodes ($\vec{f}_{\st{int}}$ in the following). \end{description} Some examples to help to understand how to use this model will be presented in the next sections. \section{Model Setup} \subsection{Setting Initial Conditions \label{sect:smm:initial_condition}} For a unique solution of the equations of motion, initial displacements and velocities for all degrees of freedom must be specified: \begin{eqnarray} \vec{u}(t=0) = \vec{u}_0\\ \dot{\vec u}(t=0) =\vec{v}_0 \end{eqnarray} The solid mechanics model can be initialized as follows: \begin{cpp} model.initFull() \end{cpp} This function initializes the internal arrays and sets them to zero. Initial displacements and velocities that are not equal to zero can be prescribed by running a loop over the total number of nodes. Here, the initial displacement in $x$-direction and the initial velocity in $y$-direction for all nodes is set to $0.1$ and $1$, respectively. \begin{cpp} auto & disp = model.getDisplacement(); auto & velo = model.getVelocity(); -for (UInt node = 0; node < mesh.getNbNodes(); ++node) { +for (Int node = 0; node < mesh.getNbNodes(); ++node) { disp(node, 0) = 0.1; velo(node, 1) = 1.; } \end{cpp} \subsection{Setting Boundary Conditions\label{sect:smm:boundary}} This section explains how to impose Dirichlet or Neumann boundary conditions. A Dirichlet boundary condition specifies the values that the displacement needs to take for every point $x$ at the boundary ($\Gamma_u$) of the problem domain (Fig.~\ref{fig:smm:boundaries}): \begin{equation} \vec{u} = \bar{\vec u} \quad \forall \vec{x}\in \Gamma_{u} \end{equation} A Neumann boundary condition imposes the value of the gradient of the solution at the boundary $\Gamma_t$ of the problem domain (Fig.~\ref{fig:smm:boundaries}): \begin{equation} \vec{t} = \mat{\sigma} \vec{n} = \bar{\vec t} \quad \forall \vec{x}\in \Gamma_{t} \end{equation} \begin{figure} \centering \def\svgwidth{0.5\columnwidth} \input{figures/problemDomain.pdf_tex} \caption{Problem domain $\Omega$ with boundary in three dimensions. The Dirchelet and the Neumann regions of the boundary are denoted with $\Gamma_u$ and $\Gamma_t$, respecitvely.\label{fig:smm:boundaries}} \label{fig:problemDomain} \end{figure} Different ways of imposing these boundary conditions exist. A basic way is to loop over nodes or elements at the boundary and apply local values. A more advanced method consists of using the notion of the boundary of the mesh. In the following both ways are presented. Starting with the basic approach, as mentioned, the Dirichlet boundary conditions can be applied by looping over the nodes and assigning the required values. Figure~\ref{fig:smm:dirichlet_bc} shows a beam with a fixed support on the left side. On the right end of the beam, a load is applied. At the fixed support, the displacement has a given value. For this example, the displacements in both the $x$ and the $y$-direction are set to zero. Implementing this displacement boundary condition is similar to the implementation of initial displacement conditions described above. However, in order to impose a displacement boundary condition for all time steps, the corresponding nodes need to be marked as boundary nodes using the function \code{blocked}. While, in order to impose a load on the right side, the nodes are not marked. The detail codes are shown as follows: \begin{cpp} auto & blocked = model.getBlockedDOFs(); const auto & pos = mesh.getNodes(); -UInt nb_nodes = mesh.getNbNodes(); +Int nb_nodes = mesh.getNbNodes(); -for (UInt node = 0; node < nb_nodes; ++node) { +for (Int node = 0; node < nb_nodes; ++node) { if(Math::are_float_equal(pos(node, _x), 0)) { blocked(node, _x) = true; // block dof in x-direction blocked(node, _y) = true; // block dof in y-direction disp(node, _x) = 0.; // fixed displacement in x-direction disp(node, _y) = 0.; // fixed displacement in y-direction } else if (Math::are_float_equal(pos(node, _y), 0)) { blocked(node, _x) = false; // unblock dof in x-direction forces(node, _x) = 10.; // force in x-direction } } \end{cpp} \begin{figure}[!htb] \centering \includegraphics[scale=0.4]{figures/dirichlet} \caption{Beam with fixed support and load.\label{fig:smm:dirichlet_bc}} \end{figure} For the more advanced approach, one needs the notion of a boundary in the mesh. Therefore, the boundary should be created before boundary condition functors can be applied. Generally the boundary can be specified from the mesh file or the geometry. For the first case, the function \code{createGroupsFromMeshData} is called. This function can read any types of mesh data which are provided in the mesh file. If the mesh file is created with Gmsh, the function takes one input strings which is either \code{tag\_0}, \code{tag\_1} or \code{physical\_names}. The first two tags are assigned by Gmsh to each element which shows the physical group that they belong to. In Gmsh, it is also possible to consider strings for different groups of elements. These elements can be separated by giving a string \code{physical\_names} to the function \code{createGroupsFromMeshData}: \begin{cpp} mesh.createGroupsFromMeshData("physical_names"). \end{cpp} Boundary conditions support can also be created from the geometry by calling \code{createBoundaryGroupFromGeometry}. This function gathers all the elements on the boundary of the geometry. To apply the required boundary conditions, the function \code{applyBC} needs to be called on a \code{SolidMechanicsModel}. This function gets a Dirichlet or Neumann functor and a string which specifies the desired boundary on which the boundary conditions is to be applied. The functors specify the type of conditions to apply. Three built-in functors for Dirichlet exist: \code{FlagOnly, FixedValue,} and \code{IncrementValue}. The functor \code{FlagOnly} is used if a point is fixed in a given direction. Therefore, the input parameter to this functor is only the fixed direction. The \code{FixedValue} functor is used when a displacement value is applied in a fixed direction. The \code{IncrementValue} applies an increment to the displacement in a given direction. The following code shows the utilization of three functors for the top, bottom and side surface of the mesh which were already defined in the Gmsh file: \begin{cpp} model.applyBC(BC::Dirichlet::FixedValue(13.0, _y), "Top"); model.applyBC(BC::Dirichlet::FlagOnly(_x), "Bottom"); model.applyBC(BC::Dirichlet::IncrementValue(13.0, _x), "Side"); \end{cpp} To apply a Neumann boundary condition, the applied traction or stress should be specified before. In case of specifying the traction on the surface, the functor \code{FromTraction} of Neumann boundary conditions is called. Otherwise, the functor \code{FromStress} should be called which gets the stress tensor as an input parameter. \begin{cpp} Vector surface_traction = {0., 0., 1.}; auto surface_stress(3, 3) = Matrix::eye(3); model.applyBC(BC::Neumann::FromTraction(surface_traction), "Bottom"); model.applyBC(BC::Neumann::FromStress(surface_stress), "Top"); \end{cpp} If the boundary conditions need to be removed during the simulation, a functor is called from the Neumann boundary condition to free those boundary conditions from the desired boundary. \begin{cpp} model.applyBC(BC::Neumann::FreeBoundary(), "Side"); \end{cpp} User specified functors can also be implemented. A full example for setting both initial and boundary conditions can be found in \shellcode{\examplesdir/boundary\_conditions.cc}. The problem solved in this example is shown in Fig.~\ref{fig:smm:bc_and_ic}. It consists of a plate that is fixed with movable supports on the left and bottom side. On the right side, a traction, which increases linearly with the number of time steps, is applied. The initial displacement and velocity in $x$-direction at all free nodes is zero and two respectively. \begin{figure}[!htb] \centering \includegraphics[scale=0.8]{figures/bc_and_ic_example} \caption{Plate on movable supports.\label{fig:smm:bc_and_ic}} \end{figure} As it is mentioned in Section \ref{sect:common:groups}, node and element groups can be used to assign the boundary conditions. A generic example is given below with a Dirichlet boundary condition. \begin{cpp} // create a node group NodeGroup & node_group = mesh.createNodeGroup("nodes_fix"); /* fill the node group with the nodes you want */ // create an element group using the existing node group mesh.createElementGroupFromNodeGroup("el_fix", "nodes_fix", spatial_dimension-1); // boundary condition can be applied using the element group name model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "el_fix"); \end{cpp} \subsection{Material Selector\label{sect:smm:materialselector}} If the user wants to assign different materials to different finite elements groups in \akantu, a material selector has to be used. By default, \akantu assigns the first valid material in the material file to all elements present in the model (regular continuum materials are assigned to the regular elements and cohesive materials are assigned to cohesive elements or element facets). To assign different materials to specific elements, mesh data information such as tag information or specified physical names can be used. \code{MeshDataMaterialSelector} class uses this information to assign different materials. With the proper physical name or tag name and index, different materials can be assigned as demonstrated in the examples below. \begin{cpp} auto mat_selector = std::make_shared>("physical_names", model); model.setMaterialSelector(mat_selector); \end{cpp} In this example the physical names specified in a GMSH geometry file will by used to match the material names in the input file. Another example would be to use the first (\code{tag\_0}) or the second (\code{tag\_1}) tag associated to each elements in the mesh: \begin{cpp} - auto mat_selector = std::make_shared>( + auto mat_selector = std::make_shared>( "tag_1", model, first_index); model.setMaterialSelector(*mat_selector); \end{cpp} where \code{first\_index} (default is 1) is the value of \code{tag\_1} that will be associated to the first material in the material input file. The following values of the tag will be associated with the following materials. There are four different material selectors pre-defined in \akantu. \code{MaterialSelector} and \code{DefaultMaterialSelector} is used to assign a material to regular elements by default. For the regular elements, as in the example above, \code{MeshDataMaterialSelector} can be used to assign different materials to different elements. Apart from the \akantu's default material selectors, users can always develop their own classes in the main code to tackle various multi-material assignment situations. % An application of \code{DefaultMaterialCohesiveSelector} and usage in % a customly generated material selector class can be seen in % \shellcode{\examplesdir/cohesive\_element/cohesive\_extrinsic\_IG\_TG/cohesive\_extrinsic\_IG\_TG.cc}. \IfFileExists{manual-cohesive_elements_insertion.tex}{\input{manual-cohesive_elements_insertion}}{} \section{Static Analysis\label{sect:smm:static}} The \code{SolidMechanicsModel} class can handle different analysis methods, the first one being presented is the static case. In this case, the equation to solve is \begin{equation} \label{eqn:smm:static} \mat{K} \vec{u} = \vec{f}_{\st{ext}} \end{equation} where $\mat{K}$ is the global stiffness matrix, $\vec{u}$ the displacement vector and $\vec{f}_{\st{ext}}$ the vector of external forces applied to the system. To solve such a problem, the static solver of the \code{SolidMechanicsModel}\index{SolidMechanicsModel} object is used. First, a model has to be created and initialized. To create the model, a mesh (which can be read from a file) is needed, as explained in Section~\ref{sect:common:mesh}. Once an instance of a \code{SolidMechanicsModel} is obtained, the easiest way to initialize it is to use the \code{initFull}\index{SolidMechanicsModel!initFull} method by giving the \code{SolidMechanicsModelOptions}. These options specify the type of analysis to be performed and whether the materials should be initialized with \code{initMaterials} or not. \begin{cpp} SolidMechanicsModel model(mesh); model.initFull(_analysis_method = _static); \end{cpp} Here, a static analysis is chosen by passing the argument \code{\_static} to the method. By default, the Boolean for no initialization of the materials is set to false, so that they are initialized during the \code{initFull}. The method \code{initFull} also initializes all appropriate vectors to zero. Once the model is created and initialized, the boundary conditions can be set as explained in Section~\ref{sect:smm:boundary}. Boundary conditions will prescribe the external forces for some free degrees of freedom $\vec{f}_{\st{ext}}$ and displacements for some others. At this point of the analysis, the function \code{solveStep}\index{SolidMechanicsModel!solveStep} can be called: \begin{cpp} model.solveStep<_scm_newton_raphson_tangent_modified, SolveConvergenceCriteria::_residual>(1e-4, 1); \end{cpp} This function is templated by the solving method and the convergence criterion and takes two arguments: the tolerance and the maximum number of iterations (100 by default), which are $\num{1e-4}$ and $1$ for this example. The modified Newton-Raphson method is chosen to solve the system. In this method, the equilibrium equation (\ref{eqn:smm:static}) is modified in order to apply a Newton-Raphson convergence algorithm: \begin{align}\label{eqn:smm:static-newton-raphson} \mat{K}^{i+1}\delta\vec{u}^{i+1} &= \vec{r} \\ &= \vec{f}_{\st{ext}} -\vec{f}_{\st{int}}\\ &= \vec{f}_{\st{ext}} - \mat{K}^{i} \vec{u}^{i}\\ \vec{u}^{i+1} &= \vec{u}^{i} + \delta\vec{u}^{i+1}~,\nonumber \end{align} where $\delta\vec{u}$ is the increment of displacement to be added from one iteration to the other, and $i$ is the Newton-Raphson iteration counter. By invoking the \code{solveStep} method in the first step, the global stiffness matrix $\mat{K}$ from Equation~(\ref{eqn:smm:static}) is automatically assembled. A Newton-Raphson iteration is subsequently started, $\mat{K}$ is updated according to the displacement computed at the previous iteration and one loops until the forces are balanced (\code{SolveConvergenceCriteria::\_residual}), \ie $||\vec{r}|| < \mbox{\code{SolveConvergenceCriteria::\_residual}}$. One can also iterate until the increment of displacement is zero (\code{SolveConvergenceCriteria::\_increment}) which also means that the equilibrium is found. For a linear elastic problem, the solution is obtained in one iteration and therefore the maximum number of iterations can be set to one. But for a non-linear case, one needs to iterate as long as the norm of the residual exceeds the tolerance threshold and therefore the maximum number of iterations has to be higher, e.g. $100$: \begin{cpp} model.solveStep<_scm_newton_raphson_tangent_modified,SolveConvergenceCriteria::_residual>(1e-4, 100) \end{cpp} At the end of the analysis, the final solution is stored in the \textbf{displacement} vector. A full example of how to solve a static problem is presented in the code \code{\examplesdir/static/static.cc}. This example is composed of a 2D plate of steel, blocked with rollers on the left and bottom sides as shown in Figure \ref{fig:smm:static}. The nodes from the right side of the sample are displaced by $0.01\%$ of the length of the plate. \begin{figure}[!htb] \centering \includegraphics[scale=1.05]{figures/static} \caption{Numerical setup\label{fig:smm:static}} \end{figure} The results of this analysis is depicted in Figure~\ref{fig:smm:implicit:static_solution}. \begin{figure}[!htb] \centering \includegraphics[width=.7\linewidth]{figures/static_analysis} \caption{Solution of the static analysis. Left: the initial condition, right: the solution (deformation magnified 50 times)} \label{fig:smm:implicit:static_solution} \end{figure} \subsection{Static implicit analysis with dynamic insertion of cohesive elements} In order to solve problems with the extrinsic cohesive method in the static implicit solution scheme, the function \code{solveStepCohesive} has to be used: \begin{cpp} model.solveStepCohesive<_scm_newton_raphson_tangent, SolveConvergenceCriteria::_increment>(1e-13, error, 25, false, 1e5, true); \end{cpp} in which the arguments are: tolerance, error, max\_iteration, load\_reduction, tol\_increase\_factor, do\_not\_factorize. This function, first applies the Newton-Raphson procedure to solve the problem. Then, it calls the method \code{checkCohesiveStress} to check if cohesive elements have to be inserted. Since the approach is implicit, only one element is added, the most stressed one (see Section \ref{extrinsic_insertion}). After insertion, the Newton-Raphson procedure is applied again to solve the same incremental loading step, with the new inserted cohesive element. The procedure loops in this way since no new cohesive elements have to be inserted. At that point, the solution is saved, and the simulation can advance to the next incremental loading step. In case the convergence is not reached, the obtained solution is not saved and the simulation return to the main file with the error given by the solution saved in the argument of the function \emph{error}. In this way, the user can intervene in the simulation in order to find anyhow convergence. A possibility is, for instance, to reduce the last incremental loading step. The variable \emph{load\_reduction} can be used to identify if the load has been already reduced or not. At the same time, with the variable \emph{tol\_increase\_factor} it is possible to increase the tolerance by a factor defined by the user in the main file, in order to accept a solution even with an error bigger than the tolerance set at the beginning. It is possible to increase the tolerance only in the phase of loading reduction, i.e., when load\_reduction = true. A not converged solution is never saved. In case the convergence is not reached even after the loading reduction procedure, the displacement field is not updated and remains the one of the last converged incremental steps. Also, cohesive elements are inserted only if convergence is reached. An example of the extrinsic cohesive method in the static implicit solution scheme is presented in \shellcode{\examplesdir/cohesive\_element/cohesive\_extrinsic\_implicit}. \section{Dynamic Methods} \label{sect:smm:Dynamic_methods} Different ways to solve the equations of motion are implemented in the solid mechanics model. The complete equations that should be solved are: \begin{equation} \label{eqn:equation-motion} \mat{M}\ddot{\vec{u}} + \mat{C}\dot{\vec{u}} + \mat{K}\vec{u} = \vec{f}_{\st{ext}}~, \end{equation} where $\mat{M}$, $\mat{C}$ and $\mat{K}$ are the mass, damping and stiffness matrices, respectively. In the previous section, it has already been discussed how to solve this equation in the static case, where $\ddot{\vec{u}} = \dot{\vec{u}} = 0$. Here the method to solve this equation in the general case will be presented. For this purpose, a time discretization has to be specified. The most common discretization method in solid mechanics is the Newmark-$\beta$ method, which is also the default in \akantu. For the Newmark-$\beta$ method, (\ref{eqn:equation-motion}) becomes a system of three equations (see \cite{curnier92a} \cite{hughes-83a} for more details): \begin{align} \mat{M} \ddot{\vec{u}}_{n+1} + \mat{C}\dot{\vec{u}}_{n+1} + \mat{K} \vec{u}_{n+1} &={\vec{f}_{\st{ext}}}_{\, n+1} \label{eqn:equation-motion-discret} \\ \vec{u}_{n+1} &=\vec{u}_{n} + \left(1 - \alpha\right) \Delta t \dot{\vec{u}}_{n} + \alpha \Delta t \dot{\vec{u}}_{n+1} + \left(\frac{1}{2} - \alpha\right) \Delta t^2 \ddot{\vec{u}}_{n} \label{eqn:finite-difference-1}\\ \dot{\vec{u}}_{n+1} &= \dot{\vec{u}}_{n} + \left(1 - \beta\right) \Delta t \ddot{\vec{u}}_{n} + \beta \Delta t \ddot{\vec{u}}_{n+1} \label{eqn:finite-difference-2} \end{align} In these new equations, $\ddot{\vec{u}}_{n}$, $\dot{\vec{u}}_{n}$ and $\vec{u}_{n}$ are the approximations of $\ddot{\vec{u}}(t_n)$, $\dot{\vec{u}}(t_n)$ and $\vec{u}(t_n)$. Equation~(\ref{eqn:equation-motion-discret}) is the equation of motion discretized in space (finite-element discretization), and equations (\ref{eqn:finite-difference-1}) and (\ref{eqn:finite-difference-2}) are discretized in both space and time (Newmark discretization). The $\alpha$ and $\beta$ parameters determine the stability and the accuracy of the algorithm. Classical values for $\alpha$ and $\beta$ are usually $\beta = 1/2$ for no numerical damping and $0 < \alpha < 1/2$. \begin{center} \begin{tabular}{cll} \toprule $\alpha$ & Method ($\beta = 1/2$) & Type\\ \midrule $0$ & central difference & explicit\\ $1/6$ & Fox-Goodwin (royal road) &implicit\\ $1/3$ & Linear acceleration &implicit\\ $1/2$ & Average acceleration (trapezoidal rule)& implicit\\ \bottomrule \end{tabular} \end{center} The solution of this system of equations, (\ref{eqn:equation-motion-discret})-(\ref{eqn:finite-difference-2}) is split into a predictor and a corrector system of equations. Moreover, in the case of a non-linear equations, an iterative algorithm such as the Newton-Raphson method is applied. The system of equations can be written as: \begin{enumerate} \item \textit{Predictor:} \begin{align} \vec{u}_{n+1}^{0} &= \vec{u}_{n} + \Delta t \dot{\vec{u}}_{n} + \frac{\Delta t^2}{2} \ddot{\vec{u}}_{n} \\ \dot{\vec{u}}_{n+1}^{0} &= \dot{\vec{u}}_{n} + \Delta t \ddot{\vec{u}}_{n} \\ \ddot{\vec{u}}_{n+1}^{0} &= \ddot{\vec{u}}_{n} \end{align} \item \textit{Solve:} \begin{align} \left(c \mat{M} + d \mat{C} + e \mat{K}_{n+1}^i\right) \vec{w} = {\vec{f}_{\st{ext}}}_{\,n+1} - {\vec{f}_{\st{int}}}_{\,n+1}^i - \mat{C} \dot{\vec{u}}_{n+1}^i - \mat{M} \ddot{\vec{u}}_{n+1}^i = \vec{r}_{n+1}^i \end{align} \item \textit{Corrector:} \begin{align} \ddot{\vec{u}}_{n+1}^{i+1} &= \ddot{\vec{u}}_{n+1}^{i} +c \vec{w} \\ \dot{\vec{u}}_{n+1}^{i+1} &= \dot{\vec{u}}_{n+1}^{i} + d\vec{w} \\ \vec{u}_{n+1}^{i+1} &= \vec{u}_{n+1}^{i} + e \vec{w} \end{align} \end{enumerate} where $i$ is the Newton-Raphson iteration counter and $c$, $d$ and $e$ are parameters depending on the method used to solve the equations \begin{center} \begin{tabular}{lcccc} \toprule & $\vec{w}$ & $e$ & $d$ & $c$\\ \midrule in acceleration &$ \delta\ddot{\vec{u}}$ & $\alpha \beta\Delta t^2$ &$\beta \Delta t$ &$1$\\ in velocity & $ \delta\dot{\vec{u}}$& $\alpha\Delta t$ & $1$ & $\frac{1}{\beta \Delta t}$\\ in displacement &$\delta\vec{u}$ & $ 1$ & $\frac{1}{\alpha \Delta t}$ & $\frac{1}{\alpha \beta \Delta t^2}$\\ \bottomrule \end{tabular} \end{center} % \note{If you want to use the implicit solver \akantu should be compiled at % least with one sparse matrix solver such as Mumps\cite{mumps}.} \subsection{Implicit Time Integration} To solve a problem with an implicit time integration scheme, first a \code{SolidMechanicsModel} object has to be created and initialized. Then the initial and boundary conditions have to be set. Everything is similar to the example in the static case (Section~\ref{sect:smm:static}), however, in this case the implicit dynamic scheme is selected at the initialization of the model. \begin{cpp} SolidMechanicsModel model(mesh); model.initFull(_analysis_method = _implicit_dynamic); /*Boundary conditions see Section ~ %\ref{sect:smm:boundary}% */ \end{cpp} Because a dynamic simulation is conducted, an integration time step $\Delta t$ has to be specified. In the case of implicit simulations, \akantu implements a trapezoidal rule by default. That is to say $\alpha = 1/2$ and $\beta = 1/2$ which is unconditionally stable. Therefore the value of the time step can be chosen arbitrarily within reason. \index{SolidMechanicsModel!setTimeStep} \begin{cpp} model.setTimeStep(time_step); \end{cpp} Since the system has to be solved for a given amount of time steps, the method \code{solveStep()}, (which has already been used in the static example in Section~\ref{sect:smm:static}), is called inside a time loop: \begin{cpp} /// time loop Real time = 0.; auto & solver = model.getNonLinearSolver(); solver.set("max_iterations", 100); solver.set("threshold", 1e-12); solver.set("convergence_type", SolveConvergenceCriteria::_solution); -for (UInt s = 1; time damage; }; \end{cpp} In order to enable to print the material parameters at any point in the user's example file using the standard output stream by typing: \begin{cpp} -for (UInt m = 0; m < model.getNbMaterials(); ++m) +for (Int m = 0; m < model.getNbMaterials(); ++m) std::cout << model.getMaterial(m) << std::endl; \end{cpp} the standard output stream operator has to be redefined. This should be done at the end of the header file: \begin{cpp} class LocalMaterialDamage : public Material { /// declare here the interace of your material }: /* ---------------------------------------------------------------------- */ /* inline functions */ /* ---------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const LocalMaterialDamage & _this) { _this.printself(stream); return stream; } \end{cpp} However, the user still needs to register the material parameters that should be printed out. The registration is done during the call of the constructor. Like all definitions the implementation of the constructor has to be written in the \code{material\_XXX.cc} file. However, the declaration has to be provided in the \code{material\_XXX.hh} file: \begin{cpp} class LocalMaterialDamage : public Material { /* -------------------------------------------------------------------- */ /* Constructors/Destructors */ /* -------------------------------------------------------------------- */ public: LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = ""); }; \end{cpp} The user can now define the implementation of the constructor in the \code{material\_XXX.cc} file: \begin{cpp} /* ---------------------------------------------------------------------- */ #include "local_material_damage.hh" #include "solid_mechanics_model.hh" namespace akantu { /* ---------------------------------------------------------------------- */ LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model, const ID & id) : Material(model, id), damage("damage", *this) { AKANTU_DEBUG_IN(); this->registerParam("E", E, 0., _pat_parsable, "Young's modulus"); this->registerParam("nu", nu, 0.5, _pat_parsable, "Poisson's ratio"); this->registerParam("lambda", lambda, _pat_readable, "First Lame coefficient"); this->registerParam("mu", mu, _pat_readable, "Second Lame coefficient"); this->registerParam("kapa", kpa, _pat_readable, "Bulk coefficient"); this->registerParam("Yd", Yd, 50., _pat_parsmod); this->registerParam("Sd", Sd, 5000., _pat_parsmod); damage.initialize(1); AKANTU_DEBUG_OUT(); } \end{cpp} During the intializer list the reference to the model and the material id are assigned and the constructor of the internal field is called. Inside the scope of the constructor the internal values have to be initialized and the parameters, that should be printed out, are registered with the function: \code{registerParam}\index{Material!registerParam}: \begin{cpp} void registerParam(name of the parameter (key in the material file), member variable, default value (optional parameter), access permissions, description); \end{cpp} The available access permissions are as follows: \begin{itemize} \item \code{\_pat\_internal}: Parameter can only be output when the material is printed. \item \code{\_pat\_writable}: User can write into the parameter. The parameter is output when the material is printed. \item \code{\_pat\_readable}: User can read the parameter. The parameter is output when the material is printed. \item \code{\_pat\_modifiable}: Parameter is writable and readable. \item \code{\_pat\_parsable}: Parameter can be parsed, \textit{i.e.} read from the input file. \item \code{\_pat\_parsmod}: Parameter is modifiable and parsable. \end{itemize} In order to implement the new constitutive law the user needs to specify how the additional material parameters, that are not defined in the input material file, should be calculated. Furthermore, it has to be defined how stresses and the stable time step should be computed for the new local material. In the case of implicit simulations, in addition, the computation of the tangent stiffness needs to be defined. Therefore, the user needs to redefine the following functions of the parent material: \begin{cpp} void initMaterial(); // for explicit and implicit simulations void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); // for implicit simulations void computeTangentStiffness(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); // for explicit and implicit simulations Real getStableTimeStep(Real h, const Element & element); \end{cpp} In the following a detailed description of these functions is provided: \begin{itemize} \item \code{initMaterial}:~ This method is called after the material file is fully read and the elements corresponding to each material are assigned. Some of the frequently used constant parameters are calculated in this method. For example, the Lam\'{e} constants of elastic materials can be considered as such parameters. \item \code{computeStress}:~ In this method, the stresses are computed based on the constitutive law as a function of the strains of the quadrature points. For example, the stresses for the elastic material are calculated based on the following formula: \begin{equation} \label{eqn:smm:constitutive_elastic} \mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu \mat{\varepsilon} \end{equation} Therefore, this method contains a loop on all quadrature points assigned to the material using the two macros:\par \code{MATERIAL\_STRESS\_QUADRATURE\_POINT\_LOOP\_BEGIN}\par \code{MATERIAL\_STRESS\_QUADRATURE\_POINT\_LOOP\_END} \begin{cpp} MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(element_type); // sigma <- f(grad_u) MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; \end{cpp} \note{The strain vector in \akantu contains the values of $\nabla \vec{u}$, i.e. it is really the \emph{displacement gradient},} \item \code{computeTangentStiffness}:~ This method is called when the tangent to the stress-strain curve is desired (see Fig \ref {fig:smm:AL:K}). For example, it is called in the implicit solver when the stiffness matrix for the regular elements is assembled based on the following formula: \begin{equation} \label{eqn:smm:constitutive_elasc} \mat{K } =\int{\mat{B^T}\mat{D(\varepsilon)}\mat{B}} \end{equation} Therefore, in this method, the \code{tangent} matrix (\mat{D}) is computed for a given strain. \note{ The \code{tangent} matrix is a $4^{th}$ order tensor which is stored as a matrix in Voigt notation.} \begin{figure}[!htb] \begin{center} \includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/tangent.pdf} \caption{Tangent to the stress-strain curve.} \label{fig:smm:AL:K} \end{center} \end{figure} \item \code{getCelerity}:~The stability criterion of the explicit integration scheme depend on the fastest wave celerity~\eqref{eqn:smm:explicit:stabletime}. This celerity depend on the material, and therefore the value of this velocity should be defined in this method for each new material. By default, the fastest wave speed is the compressive wave whose celerity can be defined in~\code{getPushWaveSpeed}. \end{itemize} Once the declaration and implementation of the new material has been completed, this material can be used in the user's example by including the header file: \begin{cpp} #include "material_XXX.hh" \end{cpp} For existing materials, as mentioned in Section~\ref{sect:smm:CL}, by default, the materials are initialized inside the method \code{initFull}. If a local material should be used instead, the initialization of the material has to be postponed until the local material is registered in the model. Therefore, the model is initialized with the boolean for skipping the material initialization equal to true: \begin{cpp} /// model initialization model.initFull(_analysis_method = _explicit_lumped_mass); \end{cpp} Once the model has been initialized, the local material needs to be registered in the model: \begin{cpp} model.registerNewCustomMaterials("name_of_local_material"); \end{cpp} Only at this point the material can be initialized: \begin{cpp} model.initMaterials(); \end{cpp} A full example for adding a new damage law can be found in \shellcode{\examplesdir/new\_material}. \subsection{Adding a New Non-Local Constitutive Law}\index{Material!create a new non-local material} In order to add a new non-local material we first have to add the local constitutive law in \akantu (see above). We can then add the non-local version of the constitutive law by adding the two files (\code{material\_XXX\_non\_local.hh} and \code{material\_XXX\_non\_local.cc}) where \code{XXX} is the name of the corresponding local material. The new law must inherit from the two classes, non-local parent class, such as the \code{MaterialNonLocal} class, and from the local version of the constitutive law, \textit{i.e.} \code{MaterialXXX}. It is therefore necessary to include the interface of those classes in the header file of your custom material and indicate the inheritance in the declaration of the class: \begin{cpp} /* ---------------------------------------------------------------------- */ #include "material_non_local.hh" // the non-local parent #include "material_XXX.hh" /* ---------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_XXX_HH__ #define __AKANTU_MATERIAL_XXX_HH__ namespace akantu { class MaterialXXXNonLocal : public MaterialXXX, public MaterialNonLocal { /// declare here the interface of your material }; \end{cpp} As members of the class we only need to add the internal fields to store the non-local quantities, which are obtained from the averaging process: \begin{cpp} /* -------------------------------------------------------------------------- */ /* Class members */ /* -------------------------------------------------------------------------- * protected: InternalField grad_u_nl; \end{cpp} The following four functions need to be implemented in the non-local material: \begin{cpp} /// initialization of the material void initMaterial(); /// loop over all element and invoke stress computation virtual void computeNonLocalStresses(GhostType ghost_type); /// compute stresses after local quantities have been averaged virtual void computeNonLocalStress(ElementType el_type, GhostType ghost_type) /// compute all local quantities void computeStress(ElementType el_type, GhostType ghost_type); \end{cpp} In the intialization of the non-local material we need to register the local quantity for the averaging process. In our example the internal field \emph{grad\_u\_nl} is the non-local counterpart of the gradient of the displacement field (\emph{grad\_u\_nl}): \begin{cpp} void MaterialXXXNonLocal::initMaterial() { MaterialXXX::initMaterial(); MaterialNonLocal::initMaterial(); /// register the non-local variable in the manager this->model->getNonLocalManager().registerNonLocalVariable(this->grad_u.getName(), this->grad_u_nl.getName(), spatial_dimension * spatial_dimension); } \end{cpp} The function to register the non-local variable takes as parameters the name of the local internal field, the name of the non-local counterpart and the number of components of the field we want to average. In the \emph{computeStress} we now need to compute all the quantities we want to average. We can then write a loop for the stress computation in the function \emph{computeNonLocalStresses} and then provide the constitutive law on each integration point in the function \emph{computeNonLocalStress}. %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/doc/manual/manual-structuralmechanicsmodel.tex b/doc/manual/manual-structuralmechanicsmodel.tex index bcf89724e..2703454f4 100644 --- a/doc/manual/manual-structuralmechanicsmodel.tex +++ b/doc/manual/manual-structuralmechanicsmodel.tex @@ -1,250 +1,250 @@ \chapter{Structural Mechanics Model} Static structural mechanics problems can be handled using the class \code{StructuralMechanicsModel}. So far, \akantu provides 2D and 3D Bernoulli beam elements \cite{frey2009}. This model is instantiated for a given \code{Mesh}, as for the \code{SolidMechanicsModel}. The model will create its own \code{FEEngine} object to compute the interpolation, gradient, integration and assembly operations. The \code{StructuralMechanicsModel} constructor is called in the following way: \begin{cpp} StructuralMechanicsModel model(mesh, spatial_dimension); \end{cpp} where \code{mesh} is a \code{Mesh} object defining the structure for which the equations of statics are to be solved, and \code{spatial\_dimension} is the dimensionality of the problem. If \code{spatial\_dimension} is omitted, the problem is assumed to have the same dimensionality as the one specified by the mesh. \note[\ 1]{Dynamic computations are not supported to date.} \note[\ 2]{Structural meshes are created and loaded as described in Section~\ref{sect:common:mesh} with \code{\_miot\_gmsh\_struct} instead of \code{\_miot\_gmsh}:} \begin{cpp} Mesh mesh; mesh.read("structural_mesh.msh", _miot_gmsh_struct); \end{cpp} This model contains at least the following \code{Arrays}: \begin{description} \item[blocked\_dofs] contains a Boolean value for each degree of freedom specifying whether that degree is blocked or not. A Dirichlet boundary condition can be prescribed by setting the \textbf{blocked\_dofs} value of a degree of freedom to \code{true}. The \textbf{displacement} is computed for all degrees of freedom for which the \textbf{blocked\_dofs} value is set to \code{false}. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept. \item[displacement\_rotation] contains the generalized displacements (\textit{i.e.} displacements and rotations) of all degrees of freedom. It can be either a computed displacement for free degrees of freedom or an imposed displacement in case of blocked ones ($\vec{u}$ in the following). \item[external\_force] contains the generalized external forces (forces and moments) applied to the nodes ($\vec{f_{\st{ext}}}$ in the following). \item[internal\_force] contains the generalized internal forces (forces and moments) applied to the nodes ($\vec{f_{\st{int}}}$ in the following). \end{description} An example to help understand how to use this model will be presented in the next section. \section{Model Setup} \label{sec:structMechMod:setup} \subsection{Initialization} The easiest way to initialize the structural mechanics model is: \begin{cpp} model.initFull(); \end{cpp} The method \code{initFull} computes the shape functions, initializes the internal vectors mentioned above and allocates the memory for the stiffness matrix, unlike the solid mechanics model, its default argument is \code{\_static}. Material properties are defined using the \code{StructuralMaterial} structure described in Table~\ref{tab:structMechMod:strucMaterial}. Such a definition could, for instance, look like \begin{cpp} StructuralMaterial mat1; mat.E=3e10; mat.I=0.0025; mat.A=0.01; \end{cpp} \begin{table}[htb] \centering \begin{tabular}{cl} \toprule Field & Description \\ \midrule \code{E} & Young's modulus \\ \code{A} & Cross section area \\ \code{I} & Second cross sectional moment of inertia (for 2D elements)\\ \code{Iy} & \code{I} around beam $y$--axis (for 3D elements)\\ \code{Iz} & \code{I} around beam $z$--axis (for 3D elements)\\ \code{GJ} & Polar moment of inertia of beam cross section (for 3D elements)\\ \bottomrule \end{tabular} \caption{Material properties for structural elements defined in the class \code{StructuralMaterial}.} \label{tab:structMechMod:strucMaterial} \end{table} Materials can be added to the model's \code{element\_material} vector using \begin{cpp} model.addMaterial(mat1); \end{cpp} They are successively numbered and then assigned to specific elements. \begin{cpp} -for (UInt i = 0; i < nb_element_mat_1; ++i) { +for (Int i = 0; i < nb_element_mat_1; ++i) { model.getElementMaterial(_bernoulli_beam_2)(i,0) = 1; } \end{cpp} \subsection{Setting Boundary Conditions}\label{sect:structMechMod:boundary} As explained before, the Dirichlet boundary conditions are applied through the array \textbf{blocked\_dofs}. Two options exist to define Neumann conditions. If a nodal force is applied, it has to be directly set in the array \textbf{force\_momentum}. For loads distributed along the beam length, the method \code{computeForcesFromFunction} integrates them into nodal forces. The method takes as input a function describing the distribution of loads along the beam and a functor \code{BoundaryFunctionType} specifing if the function is expressed in the local coordinates (\code{\_bft\_traction\_local}) or in the global system of coordinates (\code{\_bft\_traction}). \begin{cpp} static void lin_load(double * position, double * load, - Real * normal, UInt surface_id){ + Real * normal, Int surface_id){ memset(load,0,sizeof(Real)*3); load[1] = position[0]*position[0]-250; } int main(int argc, char *argv[]){ ... model.computeForcesFromFunction<_bernoulli_beam_2>(lin_load, _bft_traction_local); ...} \end{cpp} \section{Static Analysis\label{sect:structMechMod:static}} The \code{StructuralMechanicsModel} class can perform static analyses of structures. In this case, the equation to solve is the same as for the \code{SolidMechanicsModel} used for static analyses \begin{equation}\label{eqn:structMechMod:static} \mat{K} \vec{u} = \vec{f_{\st{ext}}}~, \end{equation} where $\mat{K}$ is the global stiffness matrix, $\vec{u}$ the generalized displacement vector and $\vec{f_{\st{ext}}}$ the vector of generalized external forces applied to the system. To solve such a problem, the static solver of the \code{StructuralMechanicsModel}\index{StructuralMechanicsModel} object is used. First a model has to be created and initialized. \begin{cpp} StructuralMechanicsModel model(mesh); model.initFull(); \end{cpp} \begin{itemize} \item \code{model.initFull} initializes all internal vectors to zero. \end{itemize} Once the model is created and initialized, the boundary conditions can be set as explained in Section~\ref{sect:structMechMod:boundary}. Boundary conditions will prescribe the external forces or moments for the free degrees of freedom $\vec{f_{\st{ext}}}$ and displacements or rotations for the others. To completely define the system represented by equation (\ref{eqn:structMechMod:static}), the global stiffness matrix $\mat{K}$ must be assembled. \index{StructuralMechanicsModel!assembleStiffnessMatrix} \begin{cpp} model.assembleStiffnessMatrix(); \end{cpp} The computation of the static equilibrium is performed using the same Newton-Raphson algorithm as described in Section~\ref{sect:smm:static}. \note{To date, \code{StructuralMechanicsModel} handles only constitutively and geometrically linear problems, the algorithm is therefore guaranteed to converge in two iterations.} \begin{cpp} model.updateResidual(); model.solve(); \end{cpp} \index{StructuralMechanicsModel!updateResidual} \index{StructuralMechanicsModel!solve} \begin{itemize} \item \code{model.updateResidual} assembles the internal forces and removes them from the external forces. \item \code{model.solve} solves the Equation (\ref{eqn:structMechMod:static}). The \textbf{increment} vector of the model will contain the new increment of displacements, and the \textbf{displacement\_rotation} vector is also updated to the new displacements. \end{itemize} %At the end of the analysis, the final solution is stored in the %\textbf{displacement} vector. A full example of how to solve a %structural mechanics problem is presented in the code %\shellcode{\examplesdir/structural\_mechanics/test\_structural\_mechanics\_model\_bernoulli\_beam\_2\_exemple\_1\_1.cc}. %This example is composed of a 2D beam, clamped at the left end and %supported by two rollers as shown in Figure %\ref{fig:structMechMod:exem1_1}. The problem is defined by the %applied load $q=\SI{6}{\kilo\newton\per\metre}$, moment $\bar{M} = %\SI{3.6}{\kilo\newton\metre}$, moments of inertia $I_1 = %\SI{250\,000}{\power{\centi\metre}{4}}$ and $I_2 = %\SI{128\,000}{\power{\centi\metre}{4}}$ and lengths $L_1 = %\SI{10}{\metre}$ and $L_2 = \SI{8}{\metre}$. The resulting %rotations at node two and three are $ \varphi_2 = 0.001\,167\ %\mbox{and}\ \varphi_3 = -0.000\,771.$ At the end of the analysis, the final solution is stored in the \textbf{displacement\_rotation} vector. A full example of how to solve a structural mechanics problem is presented in the code \shellcode{\examplesdir/structural\_mechanics/bernoulli\_beam\_2\_example.cc}. This example is composed of a 2D beam, clamped at the left end and supported by two rollers as shown in Figure \ref{fig:structMechMod:exem1_1}. The problem is defined by the applied load $q=\SI{6}{\kilo\newton\per\metre}$, moment $\bar{M} = \SI{3.6}{\kilo\newton\metre}$, moments of inertia $I_1 = \SI{250\,000}{\centi\metre\tothe{4}}$ and $I_2 = \SI{128\,000}{\centi\metre\tothe{4}}$ and lengths $L_1 = \SI{10}{\metre}$ and $L_2 = \SI{8}{\metre}$. The resulting rotations at node two and three are $ \varphi_2 = 0.001\,167\ \mbox{and}\ \varphi_3 = -0.000\,771.$ \begin{figure}[htb] \centering \includegraphics[scale=1.1]{figures/beam_example} \caption{2D beam example} \label{fig:structMechMod:exem1_1} \end{figure} %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh index cac962af7..ef7b989fc 100644 --- a/src/common/aka_array_tmpl.hh +++ b/src/common/aka_array_tmpl.hh @@ -1,1167 +1,1167 @@ /** * @file aka_array_tmpl.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Thu Jul 15 2010 * @date last modification: Fri Feb 26 2021 * * @brief Inline functions of the classes Array and ArrayBase * * * @section LICENSE * * Copyright (©) 2010-2021 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 . * */ /* -------------------------------------------------------------------------- */ /* Inline Functions Array */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" // NOLINT /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ //#ifndef __AKANTU_AKA_ARRAY_TMPL_HH__ //#define __AKANTU_AKA_ARRAY_TMPL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { namespace debug { struct ArrayException : public Exception {}; } // namespace debug /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(Int size, Int nb_component, const ID & id) : ArrayBase(id) { allocate(size, nb_component); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(Int size, Int nb_component, const_reference value, const ID & id) : ArrayBase(id) { allocate(size, nb_component, value); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer(const ArrayDataLayer & vect, const ID & id) : ArrayBase(vect, id) { this->data_storage = vect.data_storage; this->size_ = vect.size_; this->nb_component = vect.nb_component; this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer::ArrayDataLayer( const std::vector & vect) { this->data_storage = vect; this->size_ = vect.size(); this->nb_component = 1; this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ template ArrayDataLayer & ArrayDataLayer::operator=(const ArrayDataLayer & other) { if (this != &other) { this->data_storage = other.data_storage; this->nb_component = other.nb_component; this->size_ = other.size_; this->values = this->data_storage.data(); } return *this; } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::allocate(Int new_size, Int nb_component) { this->nb_component = nb_component; this->resize(new_size); } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::allocate(Int new_size, Int nb_component, const T & val) { this->nb_component = nb_component; this->resize(new_size, val); } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::resize(Int new_size) { this->data_storage.resize(new_size * this->nb_component); this->values = this->data_storage.data(); this->size_ = new_size; } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::resize(Int new_size, const T & value) { this->data_storage.resize(new_size * this->nb_component, value); this->values = this->data_storage.data(); this->size_ = new_size; } /* -------------------------------------------------------------------------- */ template void ArrayDataLayer::reserve(Int size, Int new_size) { if (new_size != -1) { this->data_storage.resize(new_size * this->nb_component); } this->data_storage.reserve(size * this->nb_component); this->values = this->data_storage.data(); } /* -------------------------------------------------------------------------- */ /** * append a tuple to the array with the value value for all components * @param value the new last tuple or the array will contain nb_component copies * of value */ template inline void ArrayDataLayer::push_back(const T & value) { this->data_storage.push_back(value); this->values = this->data_storage.data(); this->size_ += 1; } /* -------------------------------------------------------------------------- */ /** * append a matrix or a vector to the array * @param new_elem a reference to a Matrix or Vector */ template template inline void ArrayDataLayer::push_back( const Eigen::MatrixBase & new_elem) { AKANTU_DEBUG_ASSERT( nb_component == new_elem.size(), "The vector(" << new_elem.size() << ") as not a size compatible with the Array (nb_component=" << nb_component << ")."); for (Idx i = 0; i < new_elem.size(); ++i) { this->data_storage.push_back(new_elem.array()[i]); } this->values = this->data_storage.data(); this->size_ += 1; } /* -------------------------------------------------------------------------- */ template inline Int ArrayDataLayer::getAllocatedSize() const { return this->data_storage.capacity() / this->nb_component; } /* -------------------------------------------------------------------------- */ template inline Int ArrayDataLayer::getMemorySize() const { return this->data_storage.capacity() * sizeof(T); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class ArrayDataLayer : public ArrayBase { public: using value_type = T; using reference = value_type &; using pointer_type = value_type *; using size_type = typename ArrayBase::size_type; using const_reference = const value_type &; public: ~ArrayDataLayer() override { deallocate(); } /// Allocation of a new vector ArrayDataLayer(Int size = 0, Int nb_component = 1, const ID & id = "") : ArrayBase(id) { allocate(size, nb_component); } /// Allocation of a new vector with a default value ArrayDataLayer(Int size, Int nb_component, const_reference value, const ID & id = "") : ArrayBase(id) { allocate(size, nb_component, value); } /// Copy constructor (deep copy) ArrayDataLayer(const ArrayDataLayer & vect, const ID & id = "") : ArrayBase(vect, id) { allocate(vect.size(), vect.getNbComponent()); std::copy_n(vect.data(), this->size_ * this->nb_component, values); } /// Copy constructor (deep copy) explicit ArrayDataLayer(const std::vector & vect) { allocate(vect.size(), 1); std::copy_n(vect.data(), this->size_ * this->nb_component, values); } // copy operator inline ArrayDataLayer & operator=(const ArrayDataLayer & other) { if (this != &other) { allocate(other.size(), other.getNbComponent()); std::copy_n(other.data(), this->size_ * this->nb_component, values); } return *this; } // move constructor inline ArrayDataLayer(ArrayDataLayer && other) noexcept = default; // move assign inline ArrayDataLayer & operator=(ArrayDataLayer && other) noexcept = default; protected: // deallocate the memory virtual void deallocate() { // NOLINTNEXTLINE(cppcoreguidelines-owning-memory, // cppcoreguidelines-no-malloc) free(this->values); } // allocate the memory virtual inline void allocate(Int size, Int nb_component) { if (size != 0) { // malloc can return a non NULL pointer in case size is 0 this->values = static_cast( // NOLINT std::malloc(nb_component * size * sizeof(T))); // NOLINT } if (this->values == nullptr and size != 0) { throw std::bad_alloc(); } this->nb_component = nb_component; this->allocated_size = this->size_ = size; } // allocate and initialize the memory virtual inline void allocate(Int size, Int nb_component, const T & value) { allocate(size, nb_component); std::fill_n(values, size * nb_component, value); } public: /// append a tuple of size nb_component containing value inline void push_back(const_reference value) { resize(this->size_ + 1, value); } /// append a Vector or a Matrix template inline void push_back(const Eigen::MatrixBase & new_elem) { AKANTU_DEBUG_ASSERT( nb_component == new_elem.size(), "The vector(" << new_elem.size() << ") as not a size compatible with the Array (nb_component=" << nb_component << ")."); this->resize(this->size_ + 1); - std::copy_n(new_elem.derived().data(), new_elem.size(), - values + this->nb_component * (this->size_ - 1)); + make_view(*this, nb_component).begin()[this->size_].array() = + new_elem.array(); } /// changes the allocated size but not the size virtual void reserve(Int size, Int new_size = Int(-1)) { auto tmp_size = this->size_; if (new_size != Int(-1)) { tmp_size = new_size; } this->resize(size); this->size_ = std::min(this->size_, tmp_size); } /// change the size of the Array virtual void resize(Int size) { if (size * this->nb_component == 0) { free(values); // NOLINT: cppcoreguidelines-no-malloc values = nullptr; this->allocated_size = 0; } else { if (this->values == nullptr) { this->allocate(size, this->nb_component); return; } Int diff = size - allocated_size; Int size_to_allocate = (std::abs(diff) > AKANTU_MIN_ALLOCATION) ? size : (diff > 0) ? allocated_size + AKANTU_MIN_ALLOCATION : allocated_size; if (size_to_allocate == allocated_size) { // otherwhy the reserve + push_back might fail... this->size_ = size; return; } auto * tmp_ptr = reinterpret_cast( // NOLINT realloc(this->values, size_to_allocate * this->nb_component * sizeof(T))); if (tmp_ptr == nullptr) { throw std::bad_alloc(); } this->values = tmp_ptr; this->allocated_size = size_to_allocate; } this->size_ = size; } /// change the size of the Array and initialize the values virtual void resize(Int size, const T & val) { Int tmp_size = this->size_; this->resize(size); if (size > tmp_size) { // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) std::fill_n(values + this->nb_component * tmp_size, (size - tmp_size) * this->nb_component, val); } } /// get the amount of space allocated in bytes inline size_type getMemorySize() const final { return this->allocated_size * this->nb_component * sizeof(T); } /// Get the real size allocated in memory inline Int getAllocatedSize() const { return this->allocated_size; } /// give the address of the memory allocated for this vector [[deprecated("use data instead to be stl compatible")]] T * storage() const { return values; }; T * data() const { return values; }; protected: /// allocation type agnostic data access T * values{nullptr}; Int allocated_size{0}; }; /* -------------------------------------------------------------------------- */ template inline auto Array::operator()(Int i, Int j) -> reference { AKANTU_DEBUG_ASSERT(this->size_ > 0, "The array \"" << this->id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < this->size_) && (j < this->nb_component), "The value at position [" << i << "," << j << "] is out of range in array \"" << this->id << "\""); return this->values[i * this->nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline auto Array::operator()(Int i, Int j) const -> const_reference { AKANTU_DEBUG_ASSERT(this->size_ > 0, "The array \"" << this->id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < this->size_) && (j < this->nb_component), "The value at position [" << i << "," << j << "] is out of range in array \"" << this->id << "\""); // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) return this->values[i * this->nb_component + j]; } template inline auto Array::operator[](Int i) -> reference { AKANTU_DEBUG_ASSERT(this->size_ > 0, "The array \"" << this->id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < this->size_ * this->nb_component), "The value at position [" << i << "] is out of range in array \"" << this->id << "\""); return this->values[i]; } /* -------------------------------------------------------------------------- */ template inline auto Array::operator[](Int i) const -> const_reference { AKANTU_DEBUG_ASSERT(this->size_ > 0, "The array \"" << this->id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < this->size_ * this->nb_component), "The value at position [" << i << "] is out of range in array \"" << this->id << "\""); return this->values[i]; } /* -------------------------------------------------------------------------- */ /** * erase an element. If the erased element is not the last of the array, the * last element is moved into the hole in order to maintain contiguity. This * may invalidate existing iterators (For instance an iterator obtained by * Array::end() is no longer correct) and will change the order of the * elements. * @param i index of element to erase */ template inline void Array::erase(Idx i) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT((this->size_ > 0), "The array is empty"); AKANTU_DEBUG_ASSERT((i < this->size_), "The element at position [" << i << "] is out of range (" << i << ">=" << this->size_ << ")"); if (i != (this->size_ - 1)) { for (Idx j = 0; j < this->nb_component; ++j) { this->values[i * this->nb_component + j] = this->values[(this->size_ - 1) * this->nb_component + j]; } } this->resize(this->size_ - 1); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ template template inline auto Array::erase(const view_iterator & it) { auto && curr = it.data(); Idx pos = (curr - this->values) / this->nb_component; erase(pos); view_iterator rit = it; return --rit; } /* -------------------------------------------------------------------------- */ /** * Subtract another array entry by entry from this array in place. Both arrays * must * have the same size and nb_component. If the arrays have different shapes, * code compiled in debug mode will throw an expeption and optimised code * will behave in an unpredicted manner * @param other array to subtract from this * @return reference to modified this */ template Array & Array::operator-=(const Array & vect) { AKANTU_DEBUG_ASSERT((this->size_ == vect.size_) && (this->nb_component == vect.nb_component), "The too array don't have the same sizes"); T * a = this->values; T * b = vect.data(); for (Idx i = 0; i < this->size_ * this->nb_component; ++i) { *a -= *b; ++a; ++b; } return *this; } /* -------------------------------------------------------------------------- */ /** * Add another array entry by entry to this array in * place. Both arrays must have the same size and * nb_component. If the arrays have different shapes, code * compiled in debug mode will throw an expeption and * optimised code will behave in an unpredicted manner * @param other array to add to this * @return reference to modified this */ template Array & Array::operator+=(const Array & vect) { AKANTU_DEBUG_ASSERT((this->size_ == vect.size()) && (this->nb_component == vect.nb_component), "The too array don't have the same sizes"); T * a = this->values; T * b = vect.data(); for (Idx i = 0; i < this->size_ * this->nb_component; ++i) { *a++ += *b++; } return *this; } /* -------------------------------------------------------------------------- */ /** * Multiply all entries of this array by a scalar in place * @param alpha scalar multiplicant * @return reference to modified this */ template auto Array::operator*=(const T & alpha) -> Array & { T * a = this->values; for (Idx i = 0; i < this->size_ * this->nb_component; ++i) { *a++ *= alpha; } return *this; } /* ------------------------------------------------------------------------- */ /** * Compare this array element by element to another. * @param other array to compare to * @return true it all element are equal and arrays have * the same shape, else false */ template bool Array::operator==(const Array & other) const { bool equal = this->nb_component == other.nb_component && this->size_ == other.size_ && this->id == other.id; if (not equal) { return false; } if (this->values == other.data()) { return true; } return std::equal(this->values, this->values + this->size_ * this->nb_component, other.data()); } /* ------------------------------------------------------------------------ */ template bool Array::operator!=(const Array & other) const { return !operator==(other); } /* ------------------------------------------------------------------------ */ /** * set all tuples of the array to a given vector or matrix * @param vm Matrix or Vector to fill the array with */ template template ::value> *> inline void Array::set(const C & elem) { AKANTU_DEBUG_ASSERT( this->nb_component == elem.array().size(), "The size of the object does not match the number of components"); for (auto && v : make_view(*this, this->nb_component)) { v = elem.array(); } } /* ------------------------------------------------------------------------ */ template void Array::append(const Array & other) { AKANTU_DEBUG_ASSERT( this->nb_component == other.nb_component, "Cannot append an array with a different number of component"); Idx old_size = this->size_; this->resize(this->size_ + other.size()); T * tmp = this->values + this->nb_component * old_size; std::copy_n(other.data(), other.size() * this->nb_component, tmp); } /* ------------------------------------------------------------------------ */ /* Functions Array */ /* ------------------------------------------------------------------------ */ template Array::Array(Int size, Int nb_component, const ID & id) : parent(size, nb_component, id) {} template <> inline Array::Array(Int size, Int nb_component, const ID & id) : parent(size, nb_component, "", id) {} /* ------------------------------------------------------------------------ */ template Array::Array(Int size, Int nb_component, const_reference value, const ID & id) : parent(size, nb_component, value, id) {} /* ------------------------------------------------------------------------ */ template Array::Array(const Array & vect, const ID & id) : parent(vect, id) {} /* ------------------------------------------------------------------------ */ template auto Array::operator=(const Array & other) -> Array & { AKANTU_DEBUG_WARNING("You are copying the array " << this->id << " are you sure it is on purpose"); if (&other == this) { return *this; } parent::operator=(other); return *this; } /* ------------------------------------------------------------------------ */ template Array::Array(const std::vector & vect) : parent(vect) {} /* ------------------------------------------------------------------------ */ template Array::~Array() = default; /* ------------------------------------------------------------------------ */ /** * search elem in the array, return the position of the * first occurrence or -1 if not found * @param elem the element to look for * @return index of the first occurrence of elem or -1 if * elem is not present */ template Idx Array::find(const_reference elem) const { AKANTU_DEBUG_IN(); auto begin = this->begin(); auto end = this->end(); auto it = std::find(begin, end, elem); AKANTU_DEBUG_OUT(); return (it != end) ? it - begin : Idx(-1); } /* ------------------------------------------------------------------------ */ template template ::value> *> inline Idx Array::find(const V & elem) { AKANTU_DEBUG_ASSERT(elem.size() == this->nb_component, "Cannot find an element with a wrong size (" << elem.size() << ") != " << this->nb_component); auto && view = make_view(*this, elem.size()); auto begin = view.begin(); auto end = view.end(); auto it = std::find(begin, end, elem); return (it != end) ? it - begin : Idx(-1); } /* ------------------------------------------------------------------------ */ /** * copy the content of another array. This overwrites the * current content. * @param other Array to copy into this array. It has to * have the same nb_component as this. If compiled in * debug mode, an incorrect other will result in an * exception being thrown. Optimised code may result in * unpredicted behaviour. * @param no_sanity_check turns off all checkes */ template void Array::copy(const Array & other, bool no_sanity_check) { AKANTU_DEBUG_IN(); if (not no_sanity_check and (other.nb_component != this->nb_component)) { AKANTU_ERROR("The two arrays do not have the same " "number of components"); } this->resize((other.size_ * other.nb_component) / this->nb_component); std::copy_n(other.data(), this->size_ * this->nb_component, this->values); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ template class ArrayPrintHelper { public: template static void print_content(const Array & vect, std::ostream & stream, int indent) { std::string space(indent, AKANTU_INDENT); stream << space << " + values : {"; for (Idx i = 0; i < vect.size(); ++i) { stream << "{"; for (Idx j = 0; j < vect.getNbComponent(); ++j) { stream << vect(i, j); if (j != vect.getNbComponent() - 1) { stream << ", "; } } stream << "}"; if (i != vect.size() - 1) { stream << ", "; } } stream << "}" << std::endl; } }; template <> class ArrayPrintHelper { public: template static void print_content(__attribute__((unused)) const Array & vect, __attribute__((unused)) std::ostream & stream, __attribute__((unused)) int indent) {} }; /* ------------------------------------------------------------------------ */ template void Array::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); std::streamsize prec = stream.precision(); std::ios_base::fmtflags ff = stream.flags(); stream.setf(std::ios_base::showbase); stream.precision(2); stream << space << "Array<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; stream << space << " + id : " << this->id << std::endl; stream << space << " + size : " << this->size_ << std::endl; stream << space << " + nb_component : " << this->nb_component << std::endl; stream << space << " + allocated size : " << this->getAllocatedSize() << std::endl; stream << space << " + memory size : " << printMemorySize(this->getMemorySize()) << std::endl; if (not AKANTU_DEBUG_LEVEL_IS_TEST()) { stream << space << " + address : " << std::hex << this->values << std::dec << std::endl; } stream.precision(prec); stream.flags(ff); if (AKANTU_DEBUG_TEST(dblDump) || AKANTU_DEBUG_LEVEL_IS_TEST()) { ArrayPrintHelper::value>::print_content( *this, stream, indent); } stream << space << "]" << std::endl; } /* ------------------------------------------------------------------------ */ /* ArrayFilter */ /* ------------------------------------------------------------------------ */ template class ArrayFilter { public: /// const iterator template class const_iterator { public: Idx getCurrentIndex() { return (*this->filter_it * this->nb_item_per_elem + this->sub_element_counter); } inline const_iterator() = default; inline const_iterator(original_iterator origin_it, filter_iterator filter_it, Int nb_item_per_elem) : origin_it(std::move(origin_it)), filter_it(std::move(filter_it)), nb_item_per_elem(nb_item_per_elem), sub_element_counter(0){}; inline bool operator!=(const_iterator & other) const { return !((*this) == other); } inline bool operator==(const_iterator & other) const { return (this->origin_it == other.origin_it && this->filter_it == other.filter_it && this->sub_element_counter == other.sub_element_counter); } inline bool operator!=(const const_iterator & other) const { return !((*this) == other); } inline bool operator==(const const_iterator & other) const { return (this->origin_it == other.origin_it && this->filter_it == other.filter_it && this->sub_element_counter == other.sub_element_counter); } inline const_iterator & operator++() { ++sub_element_counter; if (sub_element_counter == nb_item_per_elem) { sub_element_counter = 0; ++filter_it; } return *this; }; inline decltype(auto) operator*() { return origin_it[nb_item_per_elem * (*filter_it) + sub_element_counter]; }; private: original_iterator origin_it; filter_iterator filter_it; /// the number of item per element Int nb_item_per_elem; /// counter for every sub element group Int sub_element_counter; }; using vector_iterator = void; // iterator>; using array_type = Array; using const_vector_iterator = const_iterator::const_scalar_iterator>; using value_type = typename array_type::value_type; private: /* ---------------------------------------------------------------------- */ /* Constructors/Destructors */ /* ---------------------------------------------------------------------- */ public: ArrayFilter(const Array & array, const Array & filter, Int nb_item_per_elem) : array(array), filter(filter), nb_item_per_elem(nb_item_per_elem){}; decltype(auto) begin_reinterpret(Int n, Int new_size) const { Int new_nb_item_per_elem = this->nb_item_per_elem; if (new_size != 0 && n != 0) new_nb_item_per_elem = this->array.getNbComponent() * this->filter.size() * this->nb_item_per_elem / (n * new_size); return const_vector_iterator(make_view(this->array, n).begin(), this->filter.begin(), new_nb_item_per_elem); }; decltype(auto) end_reinterpret(Int n, Int new_size) const { Int new_nb_item_per_elem = this->nb_item_per_elem; if (new_size != 0 && n != 0) new_nb_item_per_elem = this->array.getNbComponent() * this->filter.size() * this->nb_item_per_elem / (n * new_size); return const_vector_iterator(make_view(this->array, n).begin(), this->filter.end(), new_nb_item_per_elem); }; // vector_iterator begin_reinterpret(Int, Int) { throw; }; // vector_iterator end_reinterpret(Int, Int) { throw; }; /// return the size of the filtered array which is the filter size Int size() const { return this->filter.size() * this->nb_item_per_elem; }; /// the number of components of the filtered array Int getNbComponent() const { return this->array.getNbComponent(); }; /// tells if the container is empty bool empty() const [[gnu::warn_unused_result]] { return (size() == 0); } /* ---------------------------------------------------------------------- */ /* Class Members */ /* ---------------------------------------------------------------------- */ private: /// reference to array of data const Array & array; /// reference to the filter used to select elements const Array & filter; /// the number of item per element Int nb_item_per_elem; }; /* ------------------------------------------------------------------------ */ /* Begin/End functions implementation */ /* ------------------------------------------------------------------------ */ namespace detail { template constexpr auto take_front_impl(Tuple && t, std::index_sequence /*idxs*/) { return std::make_tuple(std::get(std::forward(t))...); } template constexpr auto take_front(Tuple && t) { return take_front_impl(std::forward(t), std::make_index_sequence{}); } template std::string to_string_all(T &&... t) { if (sizeof...(T) == 0) { return ""; } std::stringstream ss; bool noComma = true; ss << "("; (void)std::initializer_list{ (ss << (noComma ? "" : ", ") << t, noComma = false)...}; ss << ")"; return ss.str(); } template struct InstantiationHelper { template static auto instantiate(T && data, Ns... ns) { return std::make_unique(data, ns...); } }; template <> struct InstantiationHelper<0> { template static auto instantiate(T && data) { return data; } }; template decltype(auto) get_iterator(Arr && array, T * data, Ns &&... ns) { const auto is_const_arr = std::is_const>::value; using type = ViewIteratorHelper_t; using iterator = std::conditional_t, view_iterator>; static_assert(sizeof...(Ns), "You should provide a least one size"); if (array.getNbComponent() * array.size() != Int(product_all(std::forward(ns)...))) { AKANTU_CUSTOM_EXCEPTION_INFO( debug::ArrayException(), "The iterator on " << debug::demangle(typeid(Arr).name()) << to_string_all(array.size(), array.getNbComponent()) << "is not compatible with the type " << debug::demangle(typeid(type).name()) << to_string_all(ns...)); } return aka::apply([&](auto... n) { return iterator(data, n...); }, take_front(std::make_tuple(ns...))); } } // namespace detail /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ template template inline auto Array::begin(Ns &&... ns) { return detail::get_iterator(*this, this->values, std::forward(ns)..., this->size_); } template template inline auto Array::end(Ns &&... ns) { return detail::get_iterator(*this, this->values + this->nb_component * this->size_, std::forward(ns)..., this->size_); } template template inline auto Array::begin(Ns &&... ns) const { return detail::get_iterator(*this, this->values, std::forward(ns)..., this->size_); } template template inline auto Array::end(Ns &&... ns) const { return detail::get_iterator(*this, this->values + this->nb_component * this->size_, std::forward(ns)..., this->size_); } template template inline auto Array::begin_reinterpret(Ns &&... ns) { return detail::get_iterator(*this, this->values, std::forward(ns)...); } template template inline auto Array::end_reinterpret(Ns &&... ns) { return detail::get_iterator( *this, this->values + detail::product_all(std::forward(ns)...), std::forward(ns)...); } template template inline auto Array::begin_reinterpret(Ns &&... ns) const { return detail::get_iterator(*this, this->values, std::forward(ns)...); } template template inline auto Array::end_reinterpret(Ns &&... ns) const { return detail::get_iterator( *this, this->values + detail::product_all(std::forward(ns)...), std::forward(ns)...); } /* ------------------------------------------------------------------------ */ /* Views */ /* ------------------------------------------------------------------------ */ namespace detail { template class ArrayView { using tuple = std::tuple; public: using size_type = Idx; ~ArrayView() = default; constexpr ArrayView(Array && array, Ns... ns) noexcept : array(array), sizes(std::move(ns)...) {} constexpr ArrayView(const ArrayView & array_view) = default; constexpr ArrayView(ArrayView && array_view) noexcept = default; constexpr ArrayView & operator=(const ArrayView & array_view) = default; constexpr ArrayView & operator=(ArrayView && array_view) noexcept = default; auto begin() { return aka::apply( [&](auto &&... ns) { return detail::get_iterator(array.get(), array.get().data(), std::forward(ns)...); }, sizes); } auto begin() const { return aka::apply( [&](auto &&... ns) { return detail::get_iterator(array.get(), array.get().data(), std::forward(ns)...); }, sizes); } auto end() { return aka::apply( [&](auto &&... ns) { return detail::get_iterator( array.get(), array.get().data() + detail::product_all(std::forward(ns)...), std::forward(ns)...); }, sizes); } auto end() const { return aka::apply( [&](auto &&... ns) { return detail::get_iterator( array.get(), array.get().data() + detail::product_all(std::forward(ns)...), std::forward(ns)...); }, sizes); } auto cbegin() const { return this->begin(); } auto cend() const { return this->end(); } constexpr auto size() const { return std::get::value - 1>(sizes); } constexpr auto dims() const { return std::tuple_size::value - 1; } private: std::reference_wrapper> array; tuple sizes; }; /* ---------------------------------------------------------------------- */ template class ArrayView &, Ns...> { using tuple = std::tuple; public: constexpr ArrayView(const ArrayFilter & array, Ns... ns) : array(array), sizes(std::move(ns)...) {} constexpr ArrayView(const ArrayView & array_view) = default; constexpr ArrayView(ArrayView && array_view) = default; constexpr ArrayView & operator=(const ArrayView & array_view) = default; constexpr ArrayView & operator=(ArrayView && array_view) = default; auto begin() const { return aka::apply( [&](auto &&... ns) { return array.get().begin_reinterpret( std::forward(ns)...); }, sizes); } auto end() const { return aka::apply( [&](auto &&... ns) { return array.get().end_reinterpret( std::forward(ns)...); }, sizes); } auto cbegin() const { return this->begin(); } auto cend() const { return this->end(); } constexpr auto size() const { return std::get::value - 1>(sizes); } constexpr auto dims() const { return std::tuple_size::value - 1; } private: std::reference_wrapper> array; tuple sizes; }; } // namespace detail /* ------------------------------------------------------------------------ */ template decltype(auto) make_view(Array && array, const Ns... ns) { static_assert(aka::conjunction>...>::value, "Ns should be integral types"); AKANTU_DEBUG_ASSERT((detail::product_all(ns...) != 0), "You must specify non zero dimensions"); auto size = std::forward(array).size() * std::forward(array).getNbComponent() / detail::product_all(ns...); return detail::ArrayView..., std::common_type_t>( std::forward(array), std::move(ns)..., size); } } // namespace akantu //#endif /* __AKANTU_AKA_ARRAY_TMPL_HH__ */ diff --git a/src/common/aka_common.cc b/src/common/aka_common.cc index 4dca7b03b..35a14e8ec 100644 --- a/src/common/aka_common.cc +++ b/src/common/aka_common.cc @@ -1,156 +1,156 @@ /** * @file aka_common.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Wed Dec 09 2020 * * @brief Initialization of global variables * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_random_generator.hh" #include "communicator.hh" #include "cppargparse.hh" #include "parser.hh" #include "communication_tag.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace 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(); Communicator & comm = Communicator::getStaticCommunicator(argc, argv); Tag::setMaxTag(comm.getMaxTag()); debug::debugger.setParallelContext(comm.whoAmI(), comm.getNbProc()); debug::setDebugLevel(dblError); 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, (long int)(dblWarning)); static_argparser.addArgument( "--aka_print_backtrace", "Should Akantu print a backtrace in case of error", 0, cppargparse::_boolean, false, true); static_argparser.addArgument("--aka_seed", "The seed to use on prank 0", 1, cppargparse::_integer); static_argparser.parse(argc, argv, cppargparse::_remove_parsed); std::string infile = static_argparser["aka_input_file"]; if (infile.empty()) { infile = input_file; } debug::debugger.printBacktrace(static_argparser["aka_print_backtrace"]); if (not infile.empty()) { readInputFile(infile); } long int seed; if (static_argparser.has("aka_seed")) { seed = static_argparser["aka_seed"]; } else { seed = static_parser.getParameter("seed", time(nullptr), _ppsc_current_scope); } seed *= (comm.whoAmI() + 1); - RandomGenerator::seed(seed); + RandomGenerator::seed(seed); long int dbl_level = static_argparser["aka_debug_level"]; debug::setDebugLevel(DebugLevel(dbl_level)); AKANTU_DEBUG_INFO("Random seed set to " << seed); std::atexit(finalize); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void finalize() {} /* -------------------------------------------------------------------------- */ 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(ParserType::_user).first); } std::unique_ptr Communicator::static_communicator; std::ostream & operator<<(std::ostream & stream, NodeFlag flag) { using under = std::underlying_type_t; auto digits = static_cast( std::log(std::numeric_limits::max() + 1) / std::log(16)); std::ios_base::fmtflags ff; ff = stream.flags(); auto value = static_cast>(flag); stream << "0x" << std::hex << std::setw(digits) << std::setfill('0') << value; stream.flags(ff); return stream; } } // namespace akantu diff --git a/src/common/aka_extern.cc b/src/common/aka_extern.cc index d71f98cf3..e5d0913c0 100644 --- a/src/common/aka_extern.cc +++ b/src/common/aka_extern.cc @@ -1,98 +1,101 @@ /** * @file aka_extern.cc * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Tue Oct 27 2020 * * @brief initialisation of all global variables * to insure the order of creation * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_common.hh" #include "aka_math.hh" #include "aka_named_argument.hh" #include "aka_random_generator.hh" #include "communication_tag.hh" #include "cppargparse.hh" #include "parser.hh" #include "solid_mechanics_model.hh" #if defined(AKANTU_COHESIVE_ELEMENT) #include "solid_mechanics_model_cohesive.hh" #endif /* -------------------------------------------------------------------------- */ #include #include namespace akantu { /* -------------------------------------------------------------------------- */ /* error.hpp variables */ /* -------------------------------------------------------------------------- */ namespace debug { /** \todo write function to get this * values from the environment or a config file */ /// standard output for debug messages std::ostream * _akantu_debug_cout = &std::cerr; /// standard output for normal messages std::ostream & _akantu_cout = std::cout; /// parallel context used in debug messages std::string _parallel_context; Debugger debugger; } // namespace debug /* -------------------------------------------------------------------------- */ /// Paser for commandline arguments ::cppargparse::ArgumentParser static_argparser; /// Parser containing the information parsed by the input file given to initFull Parser static_parser; bool Parser::permissive_parser = false; /* -------------------------------------------------------------------------- */ Real Math::tolerance = 1e2 * std::numeric_limits::epsilon(); /* -------------------------------------------------------------------------- */ const Int _all_dimensions [[gnu::unused]] = Int(-1); /* -------------------------------------------------------------------------- */ const Array empty_filter(0, 1, "empty_filter"); /* -------------------------------------------------------------------------- */ +template <> long int RandomGenerator::_seed = 5489; +template <> std::default_random_engine RandomGenerator::generator(5489); template <> long int RandomGenerator::_seed = 5489U; template <> std::default_random_engine RandomGenerator::generator(5489U); + /* -------------------------------------------------------------------------- */ int Tag::max_tag = 0; /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/common/aka_random_generator.hh b/src/common/aka_random_generator.hh index 0312f2bd8..28e8bdb12 100644 --- a/src/common/aka_random_generator.hh +++ b/src/common/aka_random_generator.hh @@ -1,285 +1,285 @@ /** * @file aka_random_generator.hh * * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Tue Sep 29 2020 * * @brief generic random generator * * * @section LICENSE * * Copyright (©) 2014-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_AKA_RANDOM_GENERATOR_HH_ #define AKANTU_AKA_RANDOM_GENERATOR_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ /* List of available distributions */ /* -------------------------------------------------------------------------- */ // clang-format off #define AKANTU_RANDOM_DISTRIBUTION_TYPES \ ((uniform , std::uniform_real_distribution )) \ ((exponential , std::exponential_distribution )) \ ((gamma , std::gamma_distribution )) \ ((weibull , std::weibull_distribution )) \ ((extreme_value, std::extreme_value_distribution)) \ ((normal , std::normal_distribution )) \ ((lognormal , std::lognormal_distribution )) \ ((chi_squared , std::chi_squared_distribution )) \ ((cauchy , std::cauchy_distribution )) \ ((fisher_f , std::fisher_f_distribution )) \ ((student_t , std::student_t_distribution )) // clang-format on #define AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX(elem) BOOST_PP_CAT(_rdt_, elem) #define AKANTU_RANDOM_DISTRIBUTION_PREFIX(s, data, elem) \ AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX(BOOST_PP_TUPLE_ELEM(2, 0, elem)) enum RandomDistributionType { BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_RANDOM_DISTRIBUTION_PREFIX, _, AKANTU_RANDOM_DISTRIBUTION_TYPES)), _rdt_not_defined }; /* -------------------------------------------------------------------------- */ /* Generator */ /* -------------------------------------------------------------------------- */ template class RandomGenerator { /* ------------------------------------------------------------------------ */ private: static long int _seed; // NOLINT static std::default_random_engine generator; // NOLINT /* ------------------------------------------------------------------------ */ public: inline T operator()() { return generator(); } /// function to print the contain of the class void printself(std::ostream & stream, int /* indent */) const { stream << "RandGenerator [seed=" << _seed << "]"; } /* ------------------------------------------------------------------------ */ public: static void seed(long int s) { _seed = s; generator.seed(_seed); } static long int seed() { return _seed; } static constexpr T min() { return std::default_random_engine::min(); } static constexpr T max() { return std::default_random_engine::max(); } }; #if defined(__clang__) template long int RandomGenerator::_seed; // NOLINT template std::default_random_engine RandomGenerator::generator; #endif /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #undef AKANTU_RANDOM_DISTRIBUTION_PREFIX #define AKANTU_RANDOM_DISTRIBUTION_TYPE_PRINT_CASE(r, data, elem) \ case AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX( \ BOOST_PP_TUPLE_ELEM(2, 0, elem)): { \ stream << BOOST_PP_STRINGIZE(AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX( \ BOOST_PP_TUPLE_ELEM(2, 0, elem))); \ break; \ } inline std::ostream & operator<<(std::ostream & stream, RandomDistributionType type) { switch (type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_RANDOM_DISTRIBUTION_TYPE_PRINT_CASE, _, AKANTU_RANDOM_DISTRIBUTION_TYPES) default: stream << UInt(type) << " not a RandomDistributionType"; break; } return stream; } #undef AKANTU_RANDOM_DISTRIBUTION_TYPE_PRINT_CASE /* -------------------------------------------------------------------------- */ /* Some Helper */ /* -------------------------------------------------------------------------- */ template class RandomDistributionTypeHelper { enum { value = _rdt_not_defined }; }; /* -------------------------------------------------------------------------- */ #define AKANTU_RANDOM_DISTRIBUTION_TYPE_GET_TYPE(r, data, elem) \ template \ struct RandomDistributionTypeHelper > { \ enum { \ value = AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX( \ BOOST_PP_TUPLE_ELEM(2, 0, elem)) \ }; \ \ static void printself(std::ostream & stream) { \ stream << BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 0, elem)); \ } \ }; BOOST_PP_SEQ_FOR_EACH(AKANTU_RANDOM_DISTRIBUTION_TYPE_GET_TYPE, _, AKANTU_RANDOM_DISTRIBUTION_TYPES) #undef AKANTU_RANDOM_DISTRIBUTION_TYPE_GET_TYPE /* -------------------------------------------------------------------------- */ template class RandomDistribution { public: virtual ~RandomDistribution() = default; RandomDistribution() = default; RandomDistribution(const RandomDistribution & other) = default; RandomDistribution(RandomDistribution && other) noexcept = default; RandomDistribution & operator=(const RandomDistribution & other) = default; RandomDistribution & operator=(RandomDistribution && other) noexcept = default; - virtual T operator()(RandomGenerator & gen) = 0; + virtual T operator()(RandomGenerator & gen) = 0; virtual std::unique_ptr> make_unique() const = 0; virtual void printself(std::ostream & stream, int = 0) const = 0; }; template class RandomDistributionProxy : public RandomDistribution { public: explicit RandomDistributionProxy(Distribution dist) : distribution(std::move(dist)) {} - T operator()(RandomGenerator & gen) override { + T operator()(RandomGenerator & gen) override { return distribution(gen); } std::unique_ptr> make_unique() const override { return std::make_unique>( distribution); } void printself(std::ostream & stream, int /* indent */ = 0) const override { RandomDistributionTypeHelper::printself(stream); stream << " [ " << distribution << " ]"; } private: Distribution distribution; }; /* -------------------------------------------------------------------------- */ /* RandomParameter */ /* -------------------------------------------------------------------------- */ template class RandomParameter { public: template explicit RandomParameter(T base_value, Distribution dist) : base_value(base_value), type(RandomDistributionType( RandomDistributionTypeHelper::value)), distribution_proxy( std::make_unique>( std::move(dist))) {} explicit RandomParameter(T base_value) : base_value(base_value), type(RandomDistributionType( RandomDistributionTypeHelper< T, std::uniform_real_distribution>::value)), distribution_proxy( std::make_unique< RandomDistributionProxy>>( std::uniform_real_distribution(0., 0.))) {} RandomParameter(const RandomParameter & other) : base_value(other.base_value), type(other.type), distribution_proxy(other.distribution_proxy->make_unique()) {} RandomParameter & operator=(const RandomParameter & other) { distribution_proxy = other.distribution_proxy->make_unique(); base_value = other.base_value; type = other.type; return *this; } RandomParameter(RandomParameter && other) noexcept = default; RandomParameter & operator=(RandomParameter && other) noexcept = default; virtual ~RandomParameter() = default; inline void setBaseValue(const T & value) { this->base_value = value; } inline T getBaseValue() const { return this->base_value; } template