diff --git a/doc/manual/manual-feengine.tex b/doc/manual/manual-feengine.tex index 415b15d1e..ef64c823e 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, 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, 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, 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 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); +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 383b48e38..573f19639 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; 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 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{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; % 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; 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(); + 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); + 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}. +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. +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); for(UInt 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/examples/implicit/implicit_dynamic.cc b/examples/implicit/implicit_dynamic.cc index 16ed54d00..593aeca26 100644 --- a/examples/implicit/implicit_dynamic.cc +++ b/examples/implicit/implicit_dynamic.cc @@ -1,151 +1,151 @@ /** * @file implicit_dynamic.cc * * @author Nicolas Richart * * @date creation: Sun Oct 19 2014 * @date last modification: Fri Feb 28 2020 * * @brief Example of solid mechanics in implicit dynamic * * * @section LICENSE * * Copyright (©) 2015-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 "communicator.hh" #include "non_linear_solver.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ const Real bar_length = 10.; const Real bar_height = 1.; const Real bar_depth = 1.; const Real F = 5e3; const Real L = bar_length; const Real I = bar_depth * bar_height * bar_height * bar_height / 12.; const Real E = 12e7; const Real rho = 1000; const Real m = rho * bar_height * bar_depth; static Real w(UInt n) { return n * n * M_PI * M_PI / (L * L) * sqrt(E * I / m); } static Real analytical_solution(Real time) { return 2 * F * L * L * L / (pow(M_PI, 4) * E * I) * ((1. - cos(w(1) * time)) + (1. - cos(w(3) * time)) / 81. + (1. - cos(w(5) * time)) / 625.); } const UInt spatial_dimension = 2; const Real time_step = 1e-4; const Real max_time = 0.62; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("material_dynamic.dat", argc, argv); Mesh mesh(spatial_dimension); const auto & comm = Communicator::getStaticCommunicator(); Int prank = comm.whoAmI(); if (prank == 0) mesh.read("beam.msh"); mesh.distribute(); SolidMechanicsModel model(mesh); /// model initialization model.initFull(_analysis_method = _implicit_dynamic); Material & mat = model.getMaterial(0); mat.setParam("E", E); mat.setParam("rho", rho); Array & force = model.getExternalForce(); Array & displacment = model.getDisplacement(); // boundary conditions model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "blocked"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "blocked"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "roller"); - const Array & trac_nodes = + const Array & trac_nodes = mesh.getElementGroup("traction").getNodeGroup().getNodes(); bool dump_node = false; if (trac_nodes.size() > 0 && mesh.isLocalOrMasterNode(trac_nodes(0))) { force(trac_nodes(0), 1) = F; dump_node = true; } // output setup std::ofstream pos; pos.open("position.csv"); if (!pos.good()) AKANTU_ERROR("Cannot open file \"position.csv\""); pos << "id,time,position,solution" << std::endl; model.setBaseName("dynamic"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpField("external_force"); model.addDumpField("internal_force"); model.dump(); model.setTimeStep(time_step); auto & solver = model.getNonLinearSolver(); solver.set("max_iterations", 100); solver.set("threshold", 1e-12); solver.set("convergence_type", SolveConvergenceCriteria::_solution); /// time loop Real time = 0.; for (UInt s = 1; time < max_time; ++s, time += time_step) { if (prank == 0) std::cout << s << "\r" << std::flush; model.solveStep(); if (dump_node) pos << s << "," << time << "," << displacment(trac_nodes(0), 1) << "," << analytical_solution(s * time_step) << std::endl; if (s % 100 == 0) model.dump(); } std::cout << std::endl; pos.close(); finalize(); return EXIT_SUCCESS; } diff --git a/examples/io/dumper/dumpable_interface.cc b/examples/io/dumper/dumpable_interface.cc index fa5b6d2ea..41f7b92ea 100644 --- a/examples/io/dumper/dumpable_interface.cc +++ b/examples/io/dumper/dumpable_interface.cc @@ -1,191 +1,191 @@ /** * @file dumpable_interface.cc * * @author Fabian Barras * @author Nicolas Richart * * @date creation: Mon Aug 17 2015 * @date last modification: Tue Sep 29 2020 * * @brief Example usnig the dumper interface directly * * * @section LICENSE * * Copyright (©) 2015-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 "element_group.hh" #include "group_manager_inline_impl.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include "locomotive_tools.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { /* In this example, we present dumpers::Dumpable which is an interface for other classes who want to dump themselves. Several classes of Akantu inheritate from Dumpable (Model, Mesh, ...). In this example we reproduce the same tasks as example_dumper_low_level.cc using this time Dumpable interface inherted by Mesh, NodeGroup and ElementGroup. It is then advised to read first example_dumper_low_level.cc. */ initialize(argc, argv); // To start let us load the swiss train mesh and its mesh data information. UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); mesh.read("swiss_train.msh"); /* swiss_train.msh has the following physical groups that can be viewed with GMSH: "$MeshFormat 2.2 0 8 $EndMeshFormat $PhysicalNames 6 2 1 "red" 2 2 "white" 2 3 "lwheel_1" 2 4 "lwheel_2" 2 5 "rwheel_2" 2 6 "rwheel_1" $EndPhysicalNames ..." */ // Grouping nodes and elements belonging to train wheels (=four mesh data). ElementGroup & wheels_elements = mesh.createElementGroup("wheels", spatial_dimension); wheels_elements.append(mesh.getElementGroup("lwheel_1")); wheels_elements.append(mesh.getElementGroup("lwheel_2")); wheels_elements.append(mesh.getElementGroup("rwheel_1")); wheels_elements.append(mesh.getElementGroup("rwheel_2")); - const Array & lnode_1 = + const Array & lnode_1 = (mesh.getElementGroup("lwheel_1")).getNodeGroup().getNodes(); - const Array & lnode_2 = + const Array & lnode_2 = (mesh.getElementGroup("lwheel_2")).getNodeGroup().getNodes(); - const Array & rnode_1 = + const Array & rnode_1 = (mesh.getElementGroup("rwheel_1")).getNodeGroup().getNodes(); - const Array & rnode_2 = + const Array & rnode_2 = (mesh.getElementGroup("rwheel_2")).getNodeGroup().getNodes(); Array & node = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); // This time a 2D Array is created and a padding size of 3 is passed to // NodalField in order to warp train deformation on Paraview. Array displacement(nb_nodes, spatial_dimension); // Create an ElementTypeMapArray for the colour ElementTypeMapArray colour("colour"); colour.initialize(mesh, _with_nb_element = true); /* ------------------------------------------------------------------------ */ /* Creating dumpers */ /* ------------------------------------------------------------------------ */ // Create dumper for the complete mesh and register it as default dumper. auto && dumper = std::make_shared("train", "./paraview/dumpable", false); mesh.registerExternalDumper(dumper, "train", true); mesh.addDumpMesh(mesh); // The dumper for the filtered mesh can be directly taken from the // ElementGroup and then registered as "wheels_elements" dumper. auto && wheels = mesh.getGroupDumper("paraview_wheels", "wheels"); mesh.registerExternalDumper(wheels.shared_from_this(), "wheels"); mesh.setDirectoryToDumper("wheels", "./paraview/dumpable"); // Arrays and ElementTypeMapArrays can be added as external fields directly mesh.addDumpFieldExternal("displacement", displacement); ElementTypeMapArrayFilter filtered_colour( colour, wheels_elements.getElements()); auto colour_field_wheel = std::make_shared>( filtered_colour); mesh.addDumpFieldExternal("color", colour_field_wheel); mesh.addDumpFieldExternalToDumper("wheels", "displacement", displacement); mesh.addDumpFieldExternalToDumper("wheels", "colour", colour); // For some specific cases the Fields should be created, as when you want to // pad an array auto displacement_vector_field = mesh.createNodalField(&displacement, "all", 3); mesh.addDumpFieldExternal("displacement_as_paraview_vector", displacement_vector_field); mesh.addDumpFieldExternalToDumper("wheels", "displacement_as_paraview_vector", displacement_vector_field); /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ // Fill the ElementTypeMapArray colour. fillColour(mesh, colour); /// Apply displacement and wheels rotation. Real tot_displacement = 50.; Real radius = 1.; UInt nb_steps = 100; Real theta = tot_displacement / radius; Vector l_center(spatial_dimension); Vector r_center(spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { l_center(i) = node(14, i); r_center(i) = node(2, i); } for (UInt i = 0; i < nb_steps; ++i) { displacement.zero(); Real step_ratio = Real(i) / Real(nb_steps); Real angle = step_ratio * theta; applyRotation(l_center, angle, node, displacement, lnode_1); applyRotation(l_center, angle, node, displacement, lnode_2); applyRotation(r_center, angle, node, displacement, rnode_1); applyRotation(r_center, angle, node, displacement, rnode_2); for (UInt j = 0; j < nb_nodes; ++j) { displacement(j, _x) += step_ratio * tot_displacement; } /// Dump call is finally made through Dumpable interface. mesh.dump(); mesh.dump("wheels"); } finalize(); return 0; } diff --git a/examples/io/dumper/dumper_low_level.cc b/examples/io/dumper/dumper_low_level.cc index 820834668..1121dcadc 100644 --- a/examples/io/dumper/dumper_low_level.cc +++ b/examples/io/dumper/dumper_low_level.cc @@ -1,200 +1,200 @@ /** * @file dumper_low_level.cc * * @author Fabian Barras * @author Nicolas Richart * * @date creation: Mon Aug 17 2015 * @date last modification: Tue Sep 29 2020 * * @brief Example using the low level dumper interface * * * @section LICENSE * * Copyright (©) 2015-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 "element_group.hh" #include "group_manager.hh" #include "mesh.hh" #include "dumper_elemental_field.hh" #include "dumper_nodal_field.hh" #include "dumper_iohelper_paraview.hh" #include "locomotive_tools.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { /* This example aims at illustrating how to manipulate low-level methods of DumperIOHelper. The aims is to visualize a colorized moving train with Paraview */ initialize(argc, argv); // To start let us load the swiss train mesh and its mesh data information. // We aknowledge here a weel-known swiss industry for mesh donation. UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); mesh.read("swiss_train.msh"); Array & nodes = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); /* swiss_train.msh has the following physical groups that can be viewed with GMSH: "$MeshFormat 2.2 0 8 $EndMeshFormat $PhysicalNames 6 2 1 "red" 2 2 "white" 2 3 "lwheel_1" 2 4 "lwheel_2" 2 5 "rwheel_2" 2 6 "rwheel_1" $EndPhysicalNames ..." */ // Grouping nodes and elements belonging to train wheels (=four mesh data) ElementGroup & wheels_elements = mesh.createElementGroup("wheels", spatial_dimension); wheels_elements.append(mesh.getElementGroup("lwheel_1")); wheels_elements.append(mesh.getElementGroup("lwheel_2")); wheels_elements.append(mesh.getElementGroup("rwheel_1")); wheels_elements.append(mesh.getElementGroup("rwheel_2")); - const Array & lnode_1 = + const Array & lnode_1 = (mesh.getElementGroup("lwheel_1")).getNodeGroup().getNodes(); - const Array & lnode_2 = + const Array & lnode_2 = (mesh.getElementGroup("lwheel_2")).getNodeGroup().getNodes(); - const Array & rnode_1 = + const Array & rnode_1 = (mesh.getElementGroup("rwheel_1")).getNodeGroup().getNodes(); - const Array & rnode_2 = + const Array & rnode_2 = (mesh.getElementGroup("rwheel_2")).getNodeGroup().getNodes(); /* Note this Array is constructed with three components in order to warp train deformation on Paraview. A more appropriate way to do this is to set a padding in the NodalField (See example_dumpable_interface.cc.) */ Array displacement(nb_nodes, 3); // ElementalField are constructed with an ElementTypeMapArray. ElementTypeMapArray colour; colour.initialize(mesh, _with_nb_element = true); /* ------------------------------------------------------------------------ */ /* Dumper creation */ /* ------------------------------------------------------------------------ */ // Creation of two DumperParaview. One for full mesh, one for a filtered // mesh. DumperParaview dumper("train", "./paraview/dumper", false); DumperParaview wheels("wheels", "./paraview/dumper", false); // Register the full mesh dumper.registerMesh(mesh); // Register a filtered mesh limited to nodes and elements from wheels groups wheels.registerFilteredMesh(mesh, wheels_elements.getElements(), wheels_elements.getNodeGroup().getNodes()); // Generate an output file of the two mesh registered. dumper.dump(); wheels.dump(); /* At this stage no fields are attached to the two dumpers. To do so, a dumpers::Field object has to be created. Several types of dumpers::Field exist. In this example we present two of them. NodalField to describe nodal displacements of our train. ElementalField handling the color of our different part. */ // NodalField are constructed with an Array. auto displ_field = std::make_shared>(displacement); auto colour_field = std::make_shared>(colour); // Register the freshly created fields to our dumper. dumper.registerField("displacement", displ_field); dumper.registerField("colour", colour_field); // For the dumper wheels, fields have to be filtered at registration. // Filtered NodalField can be simply registered by adding an Array // listing the nodes. auto displ_field_wheel = std::make_shared>( displacement, 0, 0, &(wheels_elements.getNodeGroup().getNodes())); wheels.registerField("displacement", displ_field_wheel); // For the ElementalField, an ElementTypeMapArrayFilter has to be created. ElementTypeMapArrayFilter filtered_colour( colour, wheels_elements.getElements()); auto colour_field_wheel = std::make_shared>( filtered_colour); wheels.registerField("colour", colour_field_wheel); /* ------------------------------------------------------------------------ */ // Now that the dumpers are created and the fields are associated, let's // paint and move the train! // Fill the ElementTypeMapArray colour according to mesh data information. fillColour(mesh, colour); // Apply displacement and wheels rotation. Real tot_displacement = 50.; Real radius = 1.; UInt nb_steps = 100; Real theta = tot_displacement / radius; Vector l_center(3); Vector r_center(3); for (UInt i = 0; i < spatial_dimension; ++i) { l_center(i) = nodes(14, i); r_center(i) = nodes(2, i); } for (UInt i = 0; i < nb_steps; ++i) { displacement.zero(); Real angle = (Real)i / (Real)nb_steps * theta; applyRotation(l_center, angle, nodes, displacement, lnode_1); applyRotation(l_center, angle, nodes, displacement, lnode_2); applyRotation(r_center, angle, nodes, displacement, rnode_1); applyRotation(r_center, angle, nodes, displacement, rnode_2); for (UInt j = 0; j < nb_nodes; ++j) { displacement(j, 0) += (Real)i / (Real)nb_steps * tot_displacement; } // Output results after each moving steps for main and wheel dumpers. dumper.dump(); wheels.dump(); } finalize(); return 0; } diff --git a/examples/io/dumper/locomotive_tools.cc b/examples/io/dumper/locomotive_tools.cc index bc550fdfa..e5b26bf99 100644 --- a/examples/io/dumper/locomotive_tools.cc +++ b/examples/io/dumper/locomotive_tools.cc @@ -1,92 +1,92 @@ /** * @file locomotive_tools.cc * * @author Fabian Barras * * @date creation: Mon Aug 17 2015 * @date last modification: Fri Feb 28 2020 * * @brief Small helper code for the dumper examples * * * @section LICENSE * * Copyright (©) 2015-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 "mesh.hh" /* -------------------------------------------------------------------------- */ #include "locomotive_tools.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ void applyRotation(const Vector & center, Real angle, const Array & nodes, Array & displacement, - const Array & node_group) { + const Array & node_group) { auto nodes_it = nodes.begin(nodes.getNbComponent()); auto disp_it = displacement.begin(center.size()); - Array::const_scalar_iterator node_num_it = node_group.begin(); - Array::const_scalar_iterator node_num_end = node_group.end(); + Array::const_scalar_iterator node_num_it = node_group.begin(); + Array::const_scalar_iterator node_num_end = node_group.end(); Vector pos_rel(center.size()); for (; node_num_it != node_num_end; ++node_num_it) { const Vector pos = nodes_it[*node_num_it]; for (UInt i = 0; i < pos.size(); ++i) pos_rel(i) = pos(i); Vector dis = disp_it[*node_num_it]; pos_rel -= center; Real radius = pos_rel.norm(); if (std::abs(radius) < Math::getTolerance()) continue; Real phi_i = std::acos(pos_rel(_x) / radius); if (pos_rel(_y) < 0) phi_i *= -1; dis(_x) = std::cos(phi_i - angle) * radius - pos_rel(_x); dis(_y) = std::sin(phi_i - angle) * radius - pos_rel(_y); } } /* -------------------------------------------------------------------------- */ void fillColour(const Mesh & mesh, ElementTypeMapArray & colour) { const ElementTypeMapArray & phys_data = mesh.getData("physical_names"); const Array & txt_colour = phys_data(_triangle_3); Array & id_colour = colour(_triangle_3); for (UInt i = 0; i < txt_colour.size(); ++i) { std::string phy_name = txt_colour(i); if (phy_name == "red") id_colour(i) = 3; else if (phy_name == "white" || phy_name == "lwheel_1" || phy_name == "rwheel_1") id_colour(i) = 2; else id_colour(i) = 1; } } diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc index 87262a9f1..3b22858f2 100644 --- a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc +++ b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc @@ -1,602 +1,602 @@ /** * @file solid_mechanics_model_RVE.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 15:32:35 2016 * * @brief Implementation of SolidMechanicsModelRVE * * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_RVE.hh" #include "element_group.hh" #include "material_damage_iterative.hh" #include "node_group.hh" #include "non_linear_solver.hh" #include "non_local_manager.hh" #include "parser.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector, UInt nb_gel_pockets, UInt dim, const ID & id) : SolidMechanicsModel(mesh, dim, id), volume(0.), use_RVE_mat_selector(use_RVE_mat_selector), nb_gel_pockets(nb_gel_pockets), nb_dumps(0) { AKANTU_DEBUG_IN(); /// find the four corner nodes of the RVE findCornerNodes(); /// remove the corner nodes from the surface node groups: /// This most be done because corner nodes a not periodic mesh.getElementGroup("top").removeNode(corner_nodes(2)); mesh.getElementGroup("top").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(3)); mesh.getElementGroup("left").removeNode(corner_nodes(0)); mesh.getElementGroup("bottom").removeNode(corner_nodes(1)); mesh.getElementGroup("bottom").removeNode(corner_nodes(0)); mesh.getElementGroup("right").removeNode(corner_nodes(2)); mesh.getElementGroup("right").removeNode(corner_nodes(1)); const auto & bottom = mesh.getElementGroup("bottom").getNodeGroup(); bottom_nodes.insert(bottom.begin(), bottom.end()); const auto & left = mesh.getElementGroup("left").getNodeGroup(); left_nodes.insert(left.begin(), left.end()); // /// enforce periodicity on the displacement fluctuations // auto surface_pair_1 = std::make_pair("top", "bottom"); // auto surface_pair_2 = std::make_pair("right", "left"); // SurfacePairList surface_pairs_list; // surface_pairs_list.push_back(surface_pair_1); // surface_pairs_list.push_back(surface_pair_2); // TODO: To Nicolas correct the PBCs // this->setPBC(surface_pairs_list); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default; /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) { AKANTU_DEBUG_IN(); auto options_cp(options); options_cp.analysis_method = AnalysisMethod::_static; SolidMechanicsModel::initFullImpl(options_cp); // this->initMaterials(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); /// compute the volume of the RVE GhostType gt = _not_ghost; for (auto element_type : this->mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { Array Volume(this->mesh.getNbElement(element_type) * fem.getNbIntegrationPoints(element_type), 1, 1.); this->volume = fem.integrate(Volume, element_type); } std::cout << "The volume of the RVE is " << this->volume << std::endl; /// dumping std::stringstream base_name; base_name << this->id; this->setBaseName(base_name.str()); this->addDumpFieldVector("displacement"); this->addDumpField("stress"); this->addDumpField("grad_u"); this->addDumpField("eigen_grad_u"); this->addDumpField("blocked_dofs"); this->addDumpField("material_index"); this->addDumpField("damage"); this->addDumpField("Sc"); this->addDumpField("external_force"); this->addDumpField("equivalent_stress"); this->addDumpField("internal_force"); this->addDumpField("delta_T"); this->dump(); this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyBoundaryConditions( const Matrix & displacement_gradient) { AKANTU_DEBUG_IN(); /// get the position of the nodes const Array & pos = mesh.getNodes(); /// storage for the coordinates of a given node and the displacement that will /// be applied Vector x(spatial_dimension); Vector appl_disp(spatial_dimension); /// fix top right node UInt node = this->corner_nodes(2); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); // (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = 0.; // (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = 0.; /// apply Hx at all the other corner nodes; H: displ. gradient node = this->corner_nodes(0); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(1); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); node = this->corner_nodes(3); x(0) = pos(node, 0); x(1) = pos(node, 1); appl_disp.mul(displacement_gradient, x); (*this->blocked_dofs)(node, 0) = true; (*this->displacement)(node, 0) = appl_disp(0); (*this->blocked_dofs)(node, 1) = true; (*this->displacement)(node, 1) = appl_disp(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::applyHomogeneousTemperature( const Real & temperature) { for (UInt m = 0; m < this->getNbMaterials(); ++m) { Material & mat = this->getMaterial(m); - const ElementTypeMapArray & filter_map = mat.getElementFilter(); + const ElementTypeMapArray & filter_map = mat.getElementFilter(); // Loop over all element types for (auto && type : filter_map.elementTypes(spatial_dimension)) { - const Array & filter = filter_map(type); + const Array & filter = filter_map(type); if (filter.size() == 0) continue; Array & delta_T = mat.getArray("delta_T", type); Array::scalar_iterator delta_T_it = delta_T.begin(); Array::scalar_iterator delta_T_end = delta_T.end(); for (; delta_T_it != delta_T_end; ++delta_T_it) { *delta_T_it = temperature; } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::findCornerNodes() { AKANTU_DEBUG_IN(); // find corner nodes const auto & position = mesh.getNodes(); const auto & lower_bounds = mesh.getLowerBounds(); const auto & upper_bounds = mesh.getUpperBounds(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); corner_nodes.resize(4); corner_nodes.set(UInt(-1)); for (auto && data : enumerate(make_view(position, spatial_dimension))) { auto node = std::get<0>(data); const auto & X = std::get<1>(data); auto distance = X.distance(lower_bounds); // node 1 if (Math::are_float_equal(distance, 0)) { corner_nodes(0) = node; } // node 2 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), lower_bounds(_y))) { corner_nodes(1) = node; } // node 3 else if (Math::are_float_equal(X(_x), upper_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(2) = node; } // node 4 else if (Math::are_float_equal(X(_x), lower_bounds(_x)) && Math::are_float_equal(X(_y), upper_bounds(_y))) { corner_nodes(3) = node; } } for (UInt i = 0; i < corner_nodes.size(); ++i) { if (corner_nodes(i) == UInt(-1)) AKANTU_ERROR("The corner node " << i + 1 << " wasn't found"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::advanceASR(const Matrix & prestrain) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!"); /// apply the new eigenstrain for (auto element_type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) { Array & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = prestrain; } /// advance the damage MaterialDamageIterative<2> & mat_paste = dynamic_cast &>(*this->materials[1]); MaterialDamageIterative<2> & mat_aggregate = dynamic_cast &>(*this->materials[0]); UInt nb_damaged_elements = 0; Real max_eq_stress_aggregate = 0; Real max_eq_stress_paste = 0; auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", SolveConvergenceCriteria::_solution); do { this->solveStep(); /// compute damage max_eq_stress_aggregate = mat_aggregate.getNormMaxEquivalentStress(); max_eq_stress_paste = mat_paste.getNormMaxEquivalentStress(); nb_damaged_elements = 0; if (max_eq_stress_aggregate > max_eq_stress_paste) nb_damaged_elements = mat_aggregate.updateDamage(); else if (max_eq_stress_aggregate < max_eq_stress_paste) nb_damaged_elements = mat_paste.updateDamage(); else nb_damaged_elements = (mat_paste.updateDamage() + mat_aggregate.updateDamage()); std::cout << "the number of damaged elements is " << nb_damaged_elements << std::endl; } while (nb_damaged_elements); if (this->nb_dumps % 10 == 0) { this->dump(); } this->nb_dumps += 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine("SolidMechanicsFEEngine"); Real average = 0; GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) { if (field_type == "stress") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & stress_vec = this->materials[m]->getStress(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_stress_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_stress"); fem.integrate(stress_vec, int_stress_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) average += int_stress_vec(k, row_index * spatial_dimension + col_index); // 3 is the value // for the yy (in // 3D, the value is // 4) } } else if (field_type == "strain") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & gradu_vec = this->materials[m]->getGradU(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(gradu_vec, int_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += 0.5 * (int_gradu_vec(k, row_index * spatial_dimension + col_index) + int_gradu_vec(k, col_index * spatial_dimension + row_index)); } } else if (field_type == "eigen_grad_u") { for (UInt m = 0; m < this->materials.size(); ++m) { const auto & eigen_gradu_vec = this->materials[m]->getInternal("eigen_grad_u")(element_type); const auto & elem_filter = this->materials[m]->getElementFilter(element_type); Array int_eigen_gradu_vec(elem_filter.size(), spatial_dimension * spatial_dimension, "int_of_gradu"); fem.integrate(eigen_gradu_vec, int_eigen_gradu_vec, spatial_dimension * spatial_dimension, element_type, _not_ghost, elem_filter); for (UInt k = 0; k < elem_filter.size(); ++k) /// averaging is done only for normal components, so stress and strain /// are equal average += int_eigen_gradu_vec(k, row_index * spatial_dimension + col_index); } } else { AKANTU_ERROR("Averaging not implemented for this field!!!"); } } return average / this->volume; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeStiffness(Matrix & C_macro) { AKANTU_DEBUG_IN(); const UInt dim = 2; AKANTU_DEBUG_ASSERT(this->spatial_dimension == dim, "Is only implemented for 2D!!!"); /// apply three independent loading states to determine C /// 1. eps_el = (1;0;0) 2. eps_el = (0,1,0) 3. eps_el = (0,0,0.5) /// clear the eigenstrain Matrix zero_eigengradu(dim, dim, 0.); GhostType gt = _not_ghost; for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) { auto & prestrain_vect = const_cast &>(this->getMaterial("gel").getInternal( "eigen_grad_u")(element_type)); auto prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension); auto prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension); for (; prestrain_it != prestrain_end; ++prestrain_it) (*prestrain_it) = zero_eigengradu; } /// storage for results of 3 different loading states UInt voigt_size = VoigtHelper::size; Matrix stresses(voigt_size, voigt_size, 0.); Matrix strains(voigt_size, voigt_size, 0.); Matrix H(dim, dim, 0.); /// save the damage state before filling up cracks // ElementTypeMapReal saved_damage("saved_damage"); // saved_damage.initialize(getFEEngine(), _nb_component = 1, _default_value = // 0); // this->fillCracks(saved_damage); /// virtual test 1: H(0, 0) = 0.01; this->performVirtualTesting(H, stresses, strains, 0); /// virtual test 2: H.zero(); H(1, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 1); /// virtual test 3: H.zero(); H(0, 1) = 0.01; this->performVirtualTesting(H, stresses, strains, 2); /// drain cracks // this->drainCracks(saved_damage); /// compute effective stiffness Matrix eps_inverse(voigt_size, voigt_size); eps_inverse.inverse(strains); C_macro.mul(stresses, eps_inverse); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no) { AKANTU_DEBUG_IN(); this->applyBoundaryConditions(H); auto & solver = this->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 1e-6); solver.set("convergence_type", SolveConvergenceCriteria::_solution); this->solveStep(); /// get average stress and strain eff_stresses(0, test_no) = this->averageTensorField(0, 0, "stress"); eff_strains(0, test_no) = this->averageTensorField(0, 0, "strain"); eff_stresses(1, test_no) = this->averageTensorField(1, 1, "stress"); eff_strains(1, test_no) = this->averageTensorField(1, 1, "strain"); eff_stresses(2, test_no) = this->averageTensorField(1, 0, "stress"); eff_strains(2, test_no) = 2. * this->averageTensorField(1, 0, "strain"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::homogenizeEigenGradU( Matrix & eigen_gradu_macro) { AKANTU_DEBUG_IN(); eigen_gradu_macro(0, 0) = this->averageTensorField(0, 0, "eigen_grad_u"); eigen_gradu_macro(1, 1) = this->averageTensorField(1, 1, "eigen_grad_u"); eigen_gradu_macro(0, 1) = this->averageTensorField(0, 1, "eigen_grad_u"); eigen_gradu_macro(1, 0) = this->averageTensorField(1, 0, "eigen_grad_u"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); if (use_RVE_mat_selector) { const Vector & lowerBounds = mesh.getLowerBounds(); const Vector & upperBounds = mesh.getUpperBounds(); Real bottom = lowerBounds(1); Real top = upperBounds(1); Real box_size = std::abs(top - bottom); Real eps = box_size * 1e-6; auto tmp = std::make_shared(*this, box_size, "gel", this->nb_gel_pockets, eps); tmp->setFallback(material_selector); material_selector = tmp; } this->assignMaterialToElements(); // synchronize the element material arrays this->synchronize(SynchronizationTag::_material_id); for (auto & material : materials) { /// init internals properties const auto type = material->getID(); if (type.find("material_FE2") != std::string::npos) continue; material->initMaterial(); } this->synchronize(SynchronizationTag::_smm_init_mat); if (this->non_local_manager) { this->non_local_manager->initialize(); } // SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::fillCracks(ElementTypeMapReal & saved_damage) { const auto & mat_gel = this->getMaterial("gel"); Real E_gel = mat_gel.get("E"); Real E_homogenized = 0.; for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; Real E = mat->get("E"); auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; sav_dam = dam; for (auto q : arange(dam.size())) { E_homogenized = (E_gel - E) * dam(q) + E; dam(q) = 1. - (E_homogenized / E); } } } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelRVE::drainCracks( const ElementTypeMapReal & saved_damage) { for (auto && mat : materials) { if (mat->getName() == "gel" || mat->getName() == "FE2_mat") continue; auto & damage = mat->getInternal("damage"); for (auto && type : damage.elementTypes()) { const auto & elem_filter = mat->getElementFilter(type); auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type); auto sav_dam_it = make_view(saved_damage(type), nb_integration_point).begin(); for (auto && data : zip(elem_filter, make_view(damage(type), nb_integration_point))) { auto el = std::get<0>(data); auto & dam = std::get<1>(data); Vector sav_dam = sav_dam_it[el]; dam = sav_dam; } } } } } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh index 771229c53..516b25e59 100644 --- a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh +++ b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh @@ -1,231 +1,231 @@ /** * @file solid_mechanics_model_RVE.hh * @author Aurelia Isabel Cuba Ramos * @date Wed Jan 13 14:54:18 2016 * * @brief SMM for RVE computations in FE2 simulations * * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_RVE_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_RVE_HH_ /* -------------------------------------------------------------------------- */ #include "aka_grid_dynamic.hh" #include "solid_mechanics_model.hh" #include /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelRVE : public SolidMechanicsModel { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector = true, UInt nb_gel_pockets = 400, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model"); virtual ~SolidMechanicsModelRVE(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: void initFullImpl(const ModelOptions & option) override; /// initialize the materials void initMaterials() override; public: /// apply boundary contions based on macroscopic deformation gradient virtual void applyBoundaryConditions(const Matrix & displacement_gradient); /// apply homogeneous temperature field from the macroscale level to the RVEs virtual void applyHomogeneousTemperature(const Real & temperature); /// advance the reactions -> grow gel and apply homogenized properties void advanceASR(const Matrix & prestrain); /// compute average stress or strain in the model Real averageTensorField(UInt row_index, UInt col_index, const ID & field_type); /// compute effective stiffness of the RVE void homogenizeStiffness(Matrix & C_macro); /// compute average eigenstrain void homogenizeEigenGradU(Matrix & eigen_gradu_macro); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ inline void unpackData(CommunicationBuffer & buffer, - const Array & index, + const Array & index, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - AKANTU_GET_MACRO(CornerNodes, corner_nodes, const Array &); + AKANTU_GET_MACRO(CornerNodes, corner_nodes, const Array &); AKANTU_GET_MACRO(Volume, volume, Real); private: /// find the corner nodes void findCornerNodes(); /// perform virtual testing void performVirtualTesting(const Matrix & H, Matrix & eff_stresses, Matrix & eff_strains, const UInt test_no); void fillCracks(ElementTypeMapReal & saved_damage); void drainCracks(const ElementTypeMapReal & saved_damage); /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ /// volume of the RVE Real volume; /// corner nodes 1, 2, 3, 4 (see Leonardo's thesis, page 98) - Array corner_nodes; + Array corner_nodes; /// bottom nodes std::unordered_set bottom_nodes; /// left nodes std::unordered_set left_nodes; /// standard mat selector or user one bool use_RVE_mat_selector; /// the number of gel pockets inside the RVE UInt nb_gel_pockets; /// dump counter UInt nb_dumps; }; inline void SolidMechanicsModelRVE::unpackData(CommunicationBuffer & buffer, - const Array & index, + const Array & index, const SynchronizationTag & tag) { SolidMechanicsModel::unpackData(buffer, index, tag); // if (tag == SynchronizationTag::_smm_uv) { // auto disp_it = displacement->begin(spatial_dimension); // // for (auto node : index) { // Vector current_disp(disp_it[node]); // // // if node is at the bottom, u_bottom = u_top +u_2 -u_3 // if (bottom_nodes.count(node)) { // current_disp += Vector(disp_it[corner_nodes(1)]); // current_disp -= Vector(disp_it[corner_nodes(2)]); // } // // if node is at the left, u_left = u_right +u_4 -u_3 // else if (left_nodes.count(node)) { // current_disp += Vector(disp_it[corner_nodes(3)]); // current_disp -= Vector(disp_it[corner_nodes(2)]); // } // } // } } /* -------------------------------------------------------------------------- */ /* ASR material selector */ /* -------------------------------------------------------------------------- */ class GelMaterialSelector : public MeshDataMaterialSelector { public: GelMaterialSelector(SolidMechanicsModel & model, const Real box_size, const std::string & gel_material, const UInt nb_gel_pockets, Real /*tolerance*/ = 0.) : MeshDataMaterialSelector("physical_names", model), model(model), gel_material(gel_material), nb_gel_pockets(nb_gel_pockets), nb_placed_gel_pockets(0), box_size(box_size) { Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = model.getSpatialDimension(); Element el{_triangle_3, 0, _not_ghost}; UInt nb_element = mesh.getNbElement(el.type, el.ghost_type); Array barycenter(nb_element, spatial_dimension); for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) { el.element = std::get<0>(data); auto & bary = std::get<1>(data); mesh.getBarycenter(el, bary); } /// generate the gel pockets srand(0.); Vector center(spatial_dimension); UInt placed_gel_pockets = 0; std::set checked_baries; while (placed_gel_pockets != nb_gel_pockets) { /// get a random bary center UInt bary_id = rand() % nb_element; if (checked_baries.find(bary_id) != checked_baries.end()) continue; checked_baries.insert(bary_id); el.element = bary_id; if (MeshDataMaterialSelector::operator()(el) == 1) continue; /// element belongs to paste gel_pockets.push_back(el); placed_gel_pockets += 1; } } UInt operator()(const Element & elem) { UInt temp_index = MeshDataMaterialSelector::operator()(elem); if (temp_index == 1) return temp_index; std::vector::const_iterator iit = gel_pockets.begin(); std::vector::const_iterator eit = gel_pockets.end(); if (std::find(iit, eit, elem) != eit) { nb_placed_gel_pockets += 1; std::cout << nb_placed_gel_pockets << " gelpockets placed" << std::endl; return model.getMaterialIndex(gel_material); ; } return 0; } protected: SolidMechanicsModel & model; std::string gel_material; std::vector gel_pockets; UInt nb_gel_pockets; UInt nb_placed_gel_pockets; Real box_size; }; } // namespace akantu ///#include "material_selector_tmpl.hh" #endif /* AKANTU_SOLID_MECHANICS_MODEL_RVE_HH_ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_brittle.cc b/extra_packages/extra-materials/src/material_damage/material_brittle.cc index 70011de6f..4ca218474 100644 --- a/extra_packages/extra-materials/src/material_damage/material_brittle.cc +++ b/extra_packages/extra-materials/src/material_damage/material_brittle.cc @@ -1,123 +1,123 @@ /** * @file material_brittle.cc * * @author Aranda Ruiz Josue * @author Daniel Pino Muñoz * * * @brief Specialization of the material class for the brittle material * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_brittle.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialBrittle::MaterialBrittle(SolidMechanicsModel & model, const ID & id) : MaterialDamage(model, id), strain_rate_brittle("strain_rate_brittle", *this) { AKANTU_DEBUG_IN(); this->registerParam("S_0", S_0, 157e6, _pat_parsable | _pat_modifiable); this->registerParam("E_0", E_0, 27e3, _pat_parsable, "Strain rate threshold"); this->registerParam("A", A, 1.622e-5, _pat_parsable, "Polynome cubic constant"); this->registerParam("B", B, -1.3274, _pat_parsable, "Polynome quadratic constant"); this->registerParam("C", C, 3.6544e4, _pat_parsable, "Polynome linear constant"); this->registerParam("D", D, -181.38e6, _pat_parsable, "Polynome constant"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittle::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamage::initMaterial(); this->strain_rate_brittle.initialize(spatial_dimension * spatial_dimension); updateInternalParameters(); // this->Yd.resize(); // const Mesh & mesh = this->model.getFEEngine().getMesh(); // Mesh::type_iterator it = mesh.firstType(spatial_dimension); // Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); // for(; it != last_type; ++it) { // UInt nb_element = this->element_filter(*it).getSize(); // UInt nb_quad = this->model.getFEEngine().getNbQuadraturePoints(*it); // Array & Yd_rand_vec = Yd_rand(*it); // for(UInt e = 0; e < nb_element; ++e) { // Real rand_part = (2 * drand48()-1) * Yd_randomness * Yd; // for(UInt q = 0; q < nb_quad; ++q) // Yd_rand_vec(nb_quad*e+q,0) = Yd + rand_part; // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittle::updateInternalParameters() { MaterialDamage::updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittle::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Array & velocity = this->model.getVelocity(); Array & strain_rate_brittle = this->strain_rate_brittle(el_type, ghost_type); - Array & elem_filter = this->element_filter(el_type, ghost_type); + Array & elem_filter = this->element_filter(el_type, ghost_type); this->model.getFEEngine().gradientOnIntegrationPoints( velocity, strain_rate_brittle, spatial_dimension, el_type, ghost_type, elem_filter); Array::iterator> strain_rate_it = this->strain_rate_brittle(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Real sigma_equivalent = 0; Real fracture_stress = 0; Matrix & grad_v = *strain_rate_it; computeStressOnQuad(grad_u, grad_v, sigma, *dam, sigma_equivalent, fracture_stress); ++strain_rate_it; ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(brittle, MaterialBrittle); } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_damage/material_brittle_non_local_inline_impl.hh b/extra_packages/extra-materials/src/material_damage/material_brittle_non_local_inline_impl.hh index 8790b4350..8c164010c 100644 --- a/extra_packages/extra-materials/src/material_damage/material_brittle_non_local_inline_impl.hh +++ b/extra_packages/extra-materials/src/material_damage/material_brittle_non_local_inline_impl.hh @@ -1,116 +1,116 @@ /** * @file material_brittle_non_local_inline_impl.hh * * @author Daniel Pino Muñoz * * * @brief MaterialBrittleNonLocal inline function implementation * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ } // namespace akantu #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #include #endif namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialBrittleNonLocal::MaterialBrittleNonLocal( SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialBrittleNonLocalParent(model, id), Sigma_max("Sigma max", *this), Sigma_maxnl("Sigma_max non local", *this), Sigma_fracture("Sigma_fracture", *this) { AKANTU_DEBUG_IN(); this->is_non_local = true; this->Sigma_max.initialize(1); this->Sigma_maxnl.initialize(1); this->Sigma_fracture.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittleNonLocal::initMaterial() { AKANTU_DEBUG_IN(); this->model.getNonLocalManager().registerNonLocalVariable( this->Sigma_max.getName(), Sigma_maxnl.getName(), 1); MaterialBrittleNonLocalParent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittleNonLocal::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * Sigma_maxt = this->Sigma_max(el_type, ghost_type).storage(); Real * fracture_stress = this->Sigma_fracture(el_type, ghost_type).storage(); Array & velocity = this->model.getVelocity(); Array & strain_rate_brittle = this->strain_rate_brittle(el_type, ghost_type); - Array & elem_filter = this->element_filter(el_type, ghost_type); + Array & elem_filter = this->element_filter(el_type, ghost_type); this->model.getFEEngine().gradientOnIntegrationPoints( velocity, strain_rate_brittle, spatial_dimension, el_type, ghost_type, elem_filter); Array::iterator> strain_rate_it = this->strain_rate_brittle(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & grad_v = *strain_rate_it; MaterialBrittle::computeStressOnQuad( grad_u, grad_v, sigma, *dam, *Sigma_maxt, *fracture_stress); ++dam; ++strain_rate_it; ++Sigma_maxt; ++fracture_stress; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittleNonLocal::computeNonLocalStress( ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(type, ghost_type).storage(); Real * Sigma_maxnlt = this->Sigma_maxnl(type, ghost_type).storage(); Real * fracture_stress = this->Sigma_fracture(type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(type, ghost_type); this->computeDamageAndStressOnQuad(sigma, *dam, *Sigma_maxnlt, *fracture_stress); ++dam; ++Sigma_maxnlt; ++fracture_stress; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialBrittleNonLocal< spatial_dimension>::nonLocalVariableToNeighborhood() { this->model.getNonLocalManager().nonLocalVariableToNeighborhood( Sigma_maxnl.getName(), this->name); } diff --git a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_tmpl.hh b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_tmpl.hh index 2becb9bf1..579f68e9f 100644 --- a/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_tmpl.hh +++ b/extra_packages/extra-materials/src/material_damage/material_vreepeerlings_tmpl.hh @@ -1,106 +1,106 @@ /** * @file material_vreepeerlings_tmpl.hh * * @author Cyprien Wolff * * * @brief Specialization of the material class for the VreePeerlings material * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ template class MatParent> MaterialVreePeerlings::MaterialVreePeerlings( SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialVreePeerlingsParent(model, id), Kapa("Kapa", *this), strain_rate_vreepeerlings("strain-rate-vreepeerlings", *this), Full_dam_value_strain("fulldam-valstrain", *this), Full_dam_value_strain_rate("fulldam-valstrain-rate", *this), Number_damage("number-damage", *this), equi_strain("equi-strain", *this), equi_strain_rate("equi-strain-rate", *this) { AKANTU_DEBUG_IN(); this->registerParam("Kapaoi", Kapaoi, 0.0001, _pat_parsable); this->registerParam("Kapac", Kapac, 0.0002, _pat_parsable); this->registerParam("Arate", Arate, 0., _pat_parsable); this->registerParam("Brate", Brate, 1., _pat_parsable); this->registerParam("Crate", Brate, 1., _pat_parsable); this->registerParam("Kct", Kct, 1., _pat_parsable); this->registerParam("Kapao_randomness", Kapao_randomness, 0., _pat_parsable); this->Kapa.initialize(1); this->equi_strain.initialize(1); this->equi_strain_rate.initialize(1); this->Full_dam_value_strain.initialize(1); this->Full_dam_value_strain_rate.initialize(1); this->Number_damage.initialize(1); this->strain_rate_vreepeerlings.initialize(spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class MatParent> void MaterialVreePeerlings::initMaterial() { AKANTU_DEBUG_IN(); MaterialVreePeerlingsParent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class MatParent> void MaterialVreePeerlings::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialVreePeerlingsParent::computeStress(el_type, ghost_type); Real * dam = this->damage(el_type, ghost_type).storage(); Real * equi_straint = equi_strain(el_type, ghost_type).storage(); Real * equi_straint_rate = equi_strain_rate(el_type, ghost_type).storage(); Real * Kapaq = Kapa(el_type, ghost_type).storage(); Real * FullDam_Valstrain = Full_dam_value_strain(el_type, ghost_type).storage(); Real * FullDam_Valstrain_rate = Full_dam_value_strain_rate(el_type, ghost_type).storage(); Real * Nb_damage = Number_damage(el_type, ghost_type).storage(); Real dt = this->model.getTimeStep(); - Array & elem_filter = this->element_filter(el_type, ghost_type); + Array & elem_filter = this->element_filter(el_type, ghost_type); Array & velocity = this->model.getVelocity(); Array & strain_rate_vrplgs = this->strain_rate_vreepeerlings(el_type, ghost_type); this->model.getFEEngine().gradientOnIntegrationPoints( velocity, strain_rate_vrplgs, spatial_dimension, el_type, ghost_type, elem_filter); Array::matrix_iterator strain_rate_vrplgs_it = strain_rate_vrplgs.begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & strain_rate = *strain_rate_vrplgs_it; computeStressOnQuad(grad_u, sigma, *dam, *equi_straint, *equi_straint_rate, *Kapaq, dt, strain_rate, *FullDam_Valstrain, *FullDam_Valstrain_rate, *Nb_damage); ++dam; ++equi_straint; ++equi_straint_rate; ++Kapaq; ++strain_rate_vrplgs_it; ++FullDam_Valstrain; ++FullDam_Valstrain_rate; ++Nb_damage; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } diff --git a/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.cc b/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.cc index d6ddc6d49..bbc232c9e 100644 --- a/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.cc +++ b/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.cc @@ -1,121 +1,121 @@ /** * @file material_stiffness_proportional.cc * * @author David Simon Kammer * @author Nicolas Richart * * * @brief Special. of the material class for the caughey viscoelastic material * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_stiffness_proportional.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialStiffnessProportional::MaterialStiffnessProportional( SolidMechanicsModel & model, const ID & id) : MaterialElastic(model, id), stress_viscosity("stress_viscosity", *this), stress_elastic("stress_elastic", *this) { AKANTU_DEBUG_IN(); this->registerParam("Alpha", alpha, 0., _pat_parsable | _pat_modifiable, "Artificial viscous ratio"); this->stress_viscosity.initialize(spatial_dimension * spatial_dimension); this->stress_elastic.initialize(spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialStiffnessProportional::initMaterial() { AKANTU_DEBUG_IN(); MaterialElastic::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialStiffnessProportional::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - Array & elem_filter = this->element_filter(el_type, ghost_type); + Array & elem_filter = this->element_filter(el_type, ghost_type); Array & stress_visc = stress_viscosity(el_type, ghost_type); Array & stress_el = stress_elastic(el_type, ghost_type); MaterialElastic::computeStress(el_type, ghost_type); Array & velocity = this->model.getVelocity(); Array strain_rate(0, spatial_dimension * spatial_dimension, "strain_rate"); this->model.getFEEngine().gradientOnIntegrationPoints( velocity, strain_rate, spatial_dimension, el_type, ghost_type, elem_filter); Array::matrix_iterator strain_rate_it = strain_rate.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator stress_visc_it = stress_visc.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator stress_el_it = stress_el.begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Matrix & grad_v = *strain_rate_it; Matrix & sigma_visc = *stress_visc_it; Matrix & sigma_el = *stress_el_it; MaterialElastic::computeStressOnQuad(grad_v, sigma_visc); sigma_visc *= alpha; sigma_el.copy(sigma); sigma += sigma_visc; ++strain_rate_it; ++stress_visc_it; ++stress_el_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialStiffnessProportional::computePotentialEnergy( ElementType el_type) { AKANTU_DEBUG_IN(); Array & stress_el = stress_elastic(el_type); Array::matrix_iterator stress_el_it = stress_el.begin(spatial_dimension, spatial_dimension); Real * epot = this->potential_energy(el_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); Matrix & sigma_el = *stress_el_it; MaterialElastic::computePotentialEnergyOnQuad( grad_u, sigma_el, *epot); epot++; ++stress_el_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(ve_stiffness_prop, MaterialStiffnessProportional); } // namespace akantu diff --git a/extra_packages/extra-materials/test/test_material_damage/test_material_damage_iterative_non_local_parallel.cc b/extra_packages/extra-materials/test/test_material_damage/test_material_damage_iterative_non_local_parallel.cc index 02c42a19e..02aaf953c 100644 --- a/extra_packages/extra-materials/test/test_material_damage/test_material_damage_iterative_non_local_parallel.cc +++ b/extra_packages/extra-materials/test/test_material_damage/test_material_damage_iterative_non_local_parallel.cc @@ -1,359 +1,359 @@ /** * @file test_material_damage_iterative_non_local_parallel.cc * @author Aurelia Isabel Cuba Ramos * @date Thu Nov 26 12:20:15 2015 * * @brief test the material damage iterative non local in parallel * * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_damage_iterative.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; bool checkDisplacement(SolidMechanicsModel & model, ElementType type, std::ofstream & error_output, UInt step, bool barycenters); /* -------------------------------------------------------------------------- */ /* Main */ /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { debug::setDebugLevel(dblWarning); ElementType element_type = _triangle_3; initialize("two_materials.dat", argc, argv); const UInt spatial_dimension = 2; StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); /// read the mesh and partion it Mesh mesh(spatial_dimension); akantu::MeshPartition * partition = NULL; if (prank == 0) { mesh.read("one_circular_inclusion.msh"); /// partition the mesh partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); } /// model creation SolidMechanicsModel model(mesh); model.initParallel(partition); delete partition; /// assign the material MeshDataMaterialSelector * mat_selector; mat_selector = new MeshDataMaterialSelector("physical_names", model); model.setMaterialSelector(*mat_selector); mesh.createGroupsFromMeshData( "physical_names"); // creates groups from mesh names /// initialization of the model model.initFull(SolidMechanicsModelOptions(_static)); /// boundary conditions /// Dirichlet BC model.applyBC(BC::Dirichlet::FixedValue(0, _x), "left"); model.applyBC(BC::Dirichlet::FixedValue(0, _y), "bottom"); model.applyBC(BC::Dirichlet::FixedValue(2., _y), "top"); /// add fields that should be dumped model.setBaseName("material_damage_iterative_test"); model.addDumpFieldVector("displacement"); ; model.addDumpField("stress"); model.addDumpField("blocked_dofs"); model.addDumpField("residual"); model.addDumpField("grad_u"); model.addDumpField("damage"); model.addDumpField("partitions"); model.addDumpField("material_index"); model.addDumpField("Sc"); model.addDumpField("force"); model.addDumpField("equivalent_stress"); model.dump(); std::stringstream error_stream; error_stream << "error" << ".csv"; std::ofstream error_output; error_output.open(error_stream.str().c_str()); error_output << "# Step, Average, Max, Min" << std::endl; checkDisplacement(model, element_type, error_output, 0, true); MaterialDamageIterative & aggregate = dynamic_cast &>( model.getMaterial(0)); MaterialDamageIterative & paste = dynamic_cast &>( model.getMaterial(1)); Real error; bool converged = false; UInt nb_damaged_elements = 0; Real max_eq_stress_agg = 0; Real max_eq_stress_paste = 0; /// solve the system converged = model.solveStep<_scm_newton_raphson_tangent_modified, SolveConvergenceCriteria::_increment>(1e-12, error, 2); if (converged == false) { std::cout << "The error is: " << error << std::endl; AKANTU_DEBUG_ASSERT(converged, "Did not converge"); } if (!checkDisplacement(model, element_type, error_output, 1, false)) { finalize(); return EXIT_FAILURE; } model.dump(); /// get the maximum equivalent stress in both materials max_eq_stress_agg = aggregate.getNormMaxEquivalentStress(); max_eq_stress_paste = paste.getNormMaxEquivalentStress(); nb_damaged_elements = 0; if (max_eq_stress_agg > max_eq_stress_paste) nb_damaged_elements = aggregate.updateDamage(); else nb_damaged_elements = paste.updateDamage(); if (prank == 0 && nb_damaged_elements) std::cout << nb_damaged_elements << " elements damaged" << std::endl; /// resolve the system converged = model.solveStep<_scm_newton_raphson_tangent_modified, SolveConvergenceCriteria::_increment>(1e-12, error, 2); if (converged == false) { std::cout << "The error is: " << error << std::endl; AKANTU_DEBUG_ASSERT(converged, "Did not converge"); } if (!checkDisplacement(model, element_type, error_output, 2, false)) { finalize(); return EXIT_FAILURE; } model.dump(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ bool checkDisplacement(SolidMechanicsModel & model, ElementType type, std::ofstream & error_output, UInt step, bool barycenters) { Mesh & mesh = model.getMesh(); UInt spatial_dimension = mesh.getSpatialDimension(); - const Array & connectivity = mesh.getConnectivity(type); + const Array & connectivity = mesh.getConnectivity(type); const Array & displacement = model.getDisplacement(); UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_elem = Mesh::getNbNodesPerElement(type); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); if (psize == 1) { std::stringstream displacement_file; displacement_file << "displacement/displacement_" << std::setfill('0') << std::setw(6) << step; std::ofstream displacement_output; displacement_output.open(displacement_file.str().c_str()); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = connectivity(el, n); for (UInt dim = 0; dim < spatial_dimension; ++dim) { displacement_output << std::setprecision(15) << displacement(node, dim) << " "; } displacement_output << std::endl; } } displacement_output.close(); if (barycenters) { std::stringstream barycenter_file; barycenter_file << "displacement/barycenters"; std::ofstream barycenter_output; barycenter_output.open(barycenter_file.str().c_str()); Element element(type, 0); Vector bary(spatial_dimension); for (UInt el = 0; el < nb_element; ++el) { element.element = el; mesh.getBarycenter(element, bary); for (UInt dim = 0; dim < spatial_dimension; ++dim) { barycenter_output << std::setprecision(15) << bary(dim) << " "; } barycenter_output << std::endl; } barycenter_output.close(); } } else { if (barycenters) return true; /// read data std::stringstream displacement_file; displacement_file << "displacement/displacement_" << std::setfill('0') << std::setw(6) << step; std::ifstream displacement_input; displacement_input.open(displacement_file.str().c_str()); Array displacement_serial(0, spatial_dimension); Vector disp_tmp(spatial_dimension); while (displacement_input.good()) { for (UInt i = 0; i < spatial_dimension; ++i) displacement_input >> disp_tmp(i); displacement_serial.push_back(disp_tmp); } std::stringstream barycenter_file; barycenter_file << "displacement/barycenters"; std::ifstream barycenter_input; barycenter_input.open(barycenter_file.str().c_str()); Array barycenter_serial(0, spatial_dimension); while (barycenter_input.good()) { for (UInt dim = 0; dim < spatial_dimension; ++dim) barycenter_input >> disp_tmp(dim); barycenter_serial.push_back(disp_tmp); } Element element(type, 0); Vector bary(spatial_dimension); Array::iterator> it; Array::iterator> begin = barycenter_serial.begin(spatial_dimension); Array::iterator> end = barycenter_serial.end(spatial_dimension); Array::const_iterator> disp_it; Array::iterator> disp_serial_it; Vector difference(spatial_dimension); Array error; /// compute error for (UInt el = 0; el < nb_element; ++el) { element.element = el; mesh.getBarycenter(element, bary); /// find element for (it = begin; it != end; ++it) { UInt matched_dim = 0; while (matched_dim < spatial_dimension && Math::are_float_equal(bary(matched_dim), (*it)(matched_dim))) ++matched_dim; if (matched_dim == spatial_dimension) break; } if (it == end) { std::cout << "Element barycenter not found!" << std::endl; return false; } UInt matched_el = it - begin; disp_serial_it = displacement_serial.begin(spatial_dimension) + matched_el * nb_nodes_per_elem; for (UInt n = 0; n < nb_nodes_per_elem; ++n, ++disp_serial_it) { UInt node = connectivity(el, n); if (!mesh.isLocalOrMasterNode(node)) continue; disp_it = displacement.begin(spatial_dimension) + node; difference = *disp_it; difference -= *disp_serial_it; error.push_back(difference.norm()); } } /// compute average error Real average_error = std::accumulate(error.begin(), error.end(), 0.); comm.allReduce(&average_error, 1, _so_sum); UInt error_size = error.getSize(); comm.allReduce(&error_size, 1, _so_sum); average_error /= error_size; /// compute maximum and minimum Real max_error = *std::max_element(error.begin(), error.end()); comm.allReduce(&max_error, 1, _so_max); Real min_error = *std::min_element(error.begin(), error.end()); comm.allReduce(&min_error, 1, _so_min); /// output data if (prank == 0) { error_output << step << ", " << average_error << ", " << max_error << ", " << min_error << std::endl; } if (max_error > 1.e-9) { std::cout << "Displacement error of " << max_error << " is too big!" << std::endl; return false; } } return true; } diff --git a/extra_packages/igfem/src/fe_engine_template_tmpl_igfem.hh b/extra_packages/igfem/src/fe_engine_template_tmpl_igfem.hh index 78a9392c6..31c2db877 100644 --- a/extra_packages/igfem/src/fe_engine_template_tmpl_igfem.hh +++ b/extra_packages/igfem/src/fe_engine_template_tmpl_igfem.hh @@ -1,117 +1,117 @@ /** * @file shape_igfem_inline_impl.hh * * @author Aurelia Isabel Cuba Ramos * * @brief ShapeIGFEM inline implementation * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "integrator_gauss_igfem.hh" #include "shape_igfem.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_FE_ENGINE_TEMPLATE_TMPL_IGFEM_HH_ #define AKANTU_FE_ENGINE_TEMPLATE_TMPL_IGFEM_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ /* compatibility functions */ /* -------------------------------------------------------------------------- */ template <> inline void FEEngineTemplate:: initShapeFunctions(const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh::type_iterator it = mesh.firstType(element_dimension, ghost_type, _ek_igfem); Mesh::type_iterator end = mesh.lastType(element_dimension, ghost_type, _ek_igfem); for (; it != end; ++it) { ElementType type = *it; integrator.initIntegrator(nodes, type, ghost_type); #define INIT(_type) \ do { \ const Matrix & all_quads = \ integrator.getIntegrationPoints<_type>(ghost_type); \ const Matrix & quads_1 = integrator.getIntegrationPoints< \ ElementClassProperty<_type>::sub_element_type_1>(ghost_type); \ const Matrix & quads_2 = integrator.getIntegrationPoints< \ ElementClassProperty<_type>::sub_element_type_2>(ghost_type); \ shape_functions.initShapeFunctions(nodes, all_quads, quads_1, quads_2, \ _type, ghost_type); \ } while (0) AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(INIT); #undef INIT } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template <> inline void FEEngineTemplate:: computeIntegrationPointsCoordinates( Array & quadrature_points_coordinates, ElementType type, GhostType ghost_type, const Array & filter_elements) const { const Array & nodes_coordinates = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); /// create an array with the nodal coordinates that need to be /// interpolated. The nodal coordinates of the enriched nodes need /// to be set to zero, because they represent the enrichment of the /// position field, and the enrichments for this field are all zero! /// There is no gap in the mesh! Array igfem_nodes(nodes_coordinates.getSize(), spatial_dimension); shape_functions.extractValuesAtStandardNodes(nodes_coordinates, igfem_nodes, ghost_type); interpolateOnIntegrationPoints(igfem_nodes, quadrature_points_coordinates, spatial_dimension, type, ghost_type, filter_elements); } /* -------------------------------------------------------------------------- */ template <> inline void FEEngineTemplate:: computeIntegrationPointsCoordinates( ElementTypeMapArray & quadrature_points_coordinates, - const ElementTypeMapArray * filter_elements) const { + const ElementTypeMapArray * filter_elements) const { const Array & nodes_coordinates = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); /// create an array with the nodal coordinates that need to be /// interpolated. The nodal coordinates of the enriched nodes need /// to be set to zero, because they represent the enrichment of the /// position field, and the enrichments for this field are all zero! /// There is no gap in the mesh! Array igfem_nodes(nodes_coordinates.getSize(), spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; shape_functions.extractValuesAtStandardNodes(nodes_coordinates, igfem_nodes, ghost_type); } interpolateOnIntegrationPoints(igfem_nodes, quadrature_points_coordinates, filter_elements); } } // namespace akantu #endif /* AKANTU_FE_ENGINE_TEMPLATE_TMPL_IGFEM_HH_ */ diff --git a/extra_packages/igfem/src/material_igfem/material_igfem.cc b/extra_packages/igfem/src/material_igfem/material_igfem.cc index dceb3322f..f2d4f466a 100644 --- a/extra_packages/igfem/src/material_igfem/material_igfem.cc +++ b/extra_packages/igfem/src/material_igfem/material_igfem.cc @@ -1,367 +1,367 @@ /** * @file element_class_igfem.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * * @brief Implementation parent material for IGFEM * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ #include "material_igfem.hh" #include "aka_math.hh" namespace akantu { /* -------------------------------------------------------------------------- */ MaterialIGFEM::MaterialIGFEM(SolidMechanicsModel & model, const ID & id) : Material(model, id), nb_sub_materials(2), sub_material("sub_material", *this), name_sub_mat_1(""), name_sub_mat_2("") { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->fem = &(model.getFEEngineClass("IGFEMFEEngine")); this->model->getMesh().initElementTypeMapArray( element_filter, 1, spatial_dimension, false, _ek_igfem); this->initialize(); AKANTU_DEBUG_OUT(); }; /* -------------------------------------------------------------------------- */ MaterialIGFEM::~MaterialIGFEM() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialIGFEM::initialize() { this->gradu.setElementKind(_ek_igfem); this->stress.setElementKind(_ek_igfem); this->eigengradu.setElementKind(_ek_igfem); this->gradu.setFEEngine(*fem); this->stress.setFEEngine(*fem); this->eigengradu.setFEEngine(*fem); registerParam("name_sub_mat_1", name_sub_mat_1, std::string(), _pat_parsable | _pat_readable); registerParam("name_sub_mat_2", name_sub_mat_2, std::string(), _pat_parsable | _pat_readable); this->sub_material.initialize(1); } /* -------------------------------------------------------------------------- */ void MaterialIGFEM::computeQuadraturePointsCoordinates( ElementTypeMapArray & quadrature_points_coordinates, GhostType ghost_type) const { AKANTU_DEBUG_IN(); /// compute quadrature points position in undeformed configuration Array & nodes_coordinates = this->fem->getMesh().getNodes(); Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type, _ek_igfem); Mesh::type_iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type, _ek_igfem); for (; it != last_type; ++it) { - const Array & elem_filter = this->element_filter(*it, ghost_type); + const Array & elem_filter = this->element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element) { UInt nb_tot_quad = this->fem->getNbIntegrationPoints(*it, ghost_type) * nb_element; Array & quads = quadrature_points_coordinates(*it, ghost_type); quads.resize(nb_tot_quad); this->model->getFEEngine("IGFEMFEEngine") .interpolateOnIntegrationPoints(nodes_coordinates, quads, spatial_dimension, *it, ghost_type, elem_filter); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the gradu * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialIGFEM::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = this->fem->getMesh().firstType(spatial_dimension, ghost_type, _ek_igfem); Mesh::type_iterator last_type = this->fem->getMesh().lastType(spatial_dimension, ghost_type, _ek_igfem); for (; it != last_type; ++it) { - Array & elem_filter = element_filter(*it, ghost_type); + Array & elem_filter = element_filter(*it, ghost_type); if (elem_filter.getSize()) { Array & gradu_vect = gradu(*it, ghost_type); /// compute @f$\nabla u@f$ this->fem->gradientOnIntegrationPoints(model->getDisplacement(), gradu_vect, spatial_dimension, *it, ghost_type, elem_filter); gradu_vect -= eigengradu(*it, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(*it, ghost_type); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// extrapolate internal values void MaterialIGFEM::extrapolateInternal(const ID & id, const Element & element, const Matrix & point, Matrix & extrapolated) { if (this->isInternal(id, element.kind)) { UInt nb_element = this->element_filter(element.type, element.ghost_type).getSize(); const ID name = this->getID() + ":" + id; UInt nb_quads = this->internal_vectors_real[name]->getFEEngine().getNbIntegrationPoints( element.type, element.ghost_type); const Array & internal = this->getArray(id, element.type, element.ghost_type); UInt nb_component = internal.getNbComponent(); Array::const_matrix_iterator internal_it = internal.begin_reinterpret(nb_component, nb_quads, nb_element); Element local_element = this->convertToLocalElement(element); /// instead of really extrapolating, here the value of the first GP /// is copied into the result vector. This works only for linear /// elements /// @todo extrapolate!!!! AKANTU_DEBUG_WARNING("This is a fix, values are not truly extrapolated"); const Matrix & values = internal_it[local_element.element]; UInt index = 0; Vector tmp(nb_component); for (UInt j = 0; j < values.cols(); ++j) { tmp = values(j); if (tmp.norm() > Math::getTolerance()) { index = j; break; } } for (UInt i = 0; i < extrapolated.size(); ++i) { extrapolated(i) = values(index); } } else { Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); extrapolated = default_values; } } /* -------------------------------------------------------------------------- */ template void MaterialIGFEM::setSubMaterial(const Array & element_list, GhostType ghost_type) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template <> void MaterialIGFEM::setSubMaterial<_igfem_triangle_5>( const Array & element_list, GhostType ghost_type) { SolidMechanicsModelIGFEM * igfem_model = static_cast(this->model); Vector sub_material_index(this->nb_sub_materials); Array::const_iterator el_begin = element_list.begin(); Array::const_iterator el_end = element_list.end(); const Mesh & mesh = this->model->getMesh(); Array nodes_coordinates(mesh.getNodes(), true); Array::const_vector_iterator nodes_it = nodes_coordinates.begin(spatial_dimension); Element el; el.kind = _ek_igfem; el.type = _igfem_triangle_5; el.ghost_type = ghost_type; UInt nb_nodes_per_el = mesh.getNbNodesPerElement(el.type); UInt nb_parent_nodes = IGFEMHelper::getNbParentNodes(el.type); Vector is_inside(nb_parent_nodes); - const Array & connectivity = mesh.getConnectivity(el.type, ghost_type); - Array::const_vector_iterator connec_it = + const Array & connectivity = mesh.getConnectivity(el.type, ghost_type); + Array::const_vector_iterator connec_it = connectivity.begin(nb_nodes_per_el); /// get the number of quadrature points for the two sub-elements UInt quads_1 = IGFEMHelper::getNbQuadraturePoints(el.type, 0); UInt quads_2 = IGFEMHelper::getNbQuadraturePoints(el.type, 1); UInt nb_total_quads = quads_1 + quads_2; UInt * sub_mat_ptr = this->sub_material(el.type, ghost_type).storage(); /// loop all elements for the given type - const Array & filter = this->element_filter(el.type, ghost_type); + const Array & filter = this->element_filter(el.type, ghost_type); UInt nb_elements = filter.getSize(); for (UInt e = 0; e < nb_elements; ++e, ++connec_it) { el.element = filter(e); if (std::find(el_begin, el_end, el) == el_end) { sub_mat_ptr += nb_total_quads; continue; } for (UInt i = 0; i < nb_parent_nodes; ++i) { Vector node = nodes_it[(*connec_it)(i)]; is_inside(i) = igfem_model->isInside(node, this->name_sub_mat_1); } UInt orientation = IGFEMHelper::getElementOrientation(el.type, is_inside); switch (orientation) { case 0: { sub_material_index(0) = 0; sub_material_index(1) = 1; break; } case 1: { sub_material_index(0) = 1; sub_material_index(1) = 0; break; } case 2: { sub_material_index(0) = 0; sub_material_index(1) = 0; break; } case 3: { sub_material_index(0) = 1; sub_material_index(0) = 1; break; } } for (UInt q = 0; q < quads_1; ++q, ++sub_mat_ptr) { UInt index = sub_material_index(0); *sub_mat_ptr = index; } for (UInt q = 0; q < quads_2; ++q, ++sub_mat_ptr) { UInt index = sub_material_index(1); *sub_mat_ptr = index; } } } /* -------------------------------------------------------------------------- */ template <> void MaterialIGFEM::setSubMaterial<_igfem_triangle_4>( const Array & element_list, GhostType ghost_type) { SolidMechanicsModelIGFEM * igfem_model = static_cast(this->model); Vector sub_material_index(this->nb_sub_materials); Array::const_iterator el_begin = element_list.begin(); Array::const_iterator el_end = element_list.end(); const Mesh & mesh = this->model->getMesh(); Element el; el.kind = _ek_igfem; el.ghost_type = ghost_type; el.type = _igfem_triangle_4; UInt nb_nodes_per_el = mesh.getNbNodesPerElement(el.type); Vector barycenter(spatial_dimension); - const Array & connectivity = mesh.getConnectivity(el.type, ghost_type); - Array::const_vector_iterator connec_it = + const Array & connectivity = mesh.getConnectivity(el.type, ghost_type); + Array::const_vector_iterator connec_it = connectivity.begin(nb_nodes_per_el); /// get the number of quadrature points for the two sub-elements UInt quads_1 = IGFEMHelper::getNbQuadraturePoints(el.type, 0); UInt quads_2 = IGFEMHelper::getNbQuadraturePoints(el.type, 1); UInt nb_total_quads = quads_1 + quads_2; UInt * sub_mat_ptr = this->sub_material(el.type, ghost_type).storage(); /// loop all elements for the given type - const Array & filter = this->element_filter(el.type, ghost_type); + const Array & filter = this->element_filter(el.type, ghost_type); UInt nb_elements = filter.getSize(); for (UInt e = 0; e < nb_elements; ++e, ++connec_it) { el.element = filter(e); if (std::find(el_begin, el_end, el) == el_end) { sub_mat_ptr += nb_total_quads; continue; } for (UInt s = 0; s < this->nb_sub_materials; ++s) { igfem_model->getSubElementBarycenter(el.element, s, el.type, barycenter, ghost_type); sub_material_index(s) = 1 - igfem_model->isInside(barycenter, this->name_sub_mat_1); } for (UInt q = 0; q < quads_1; ++q, ++sub_mat_ptr) { UInt index = sub_material_index(0); *sub_mat_ptr = index; } for (UInt q = 0; q < quads_2; ++q, ++sub_mat_ptr) { UInt index = sub_material_index(1); *sub_mat_ptr = index; } } } /* -------------------------------------------------------------------------- */ void MaterialIGFEM::applyEigenGradU( const Matrix & prescribed_eigen_grad_u, const ID & id, const GhostType ghost_type) { std::map::const_iterator sub_mat_it = this->sub_material_names.begin(); for (; sub_mat_it != sub_material_names.end(); ++sub_mat_it) { if (sub_mat_it->second == id) { UInt sub_element_index = sub_mat_it->first; - ElementTypeMapArray::type_iterator it = + ElementTypeMapArray::type_iterator it = this->element_filter.firstType(_all_dimensions, ghost_type, _ek_not_defined); - ElementTypeMapArray::type_iterator end = + ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, ghost_type, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if (!element_filter(type, ghost_type).getSize()) continue; Array::matrix_iterator eigen_it = this->eigengradu(type, ghost_type) .begin(spatial_dimension, spatial_dimension); Array::matrix_iterator eigen_end = this->eigengradu(type, ghost_type) .end(spatial_dimension, spatial_dimension); UInt * sub_mat_ptr = this->sub_material(type, ghost_type).storage(); for (; eigen_it != eigen_end; ++eigen_it, ++sub_mat_ptr) { if (*sub_mat_ptr == sub_element_index) { Matrix & current_eigengradu = *eigen_it; current_eigengradu = prescribed_eigen_grad_u; } } } } } } } // namespace akantu diff --git a/extra_packages/igfem/src/material_igfem/material_igfem.hh b/extra_packages/igfem/src/material_igfem/material_igfem.hh index 1d708d2ec..fe3d02c04 100644 --- a/extra_packages/igfem/src/material_igfem/material_igfem.hh +++ b/extra_packages/igfem/src/material_igfem/material_igfem.hh @@ -1,134 +1,134 @@ /** * @file element_class_igfem.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * * @brief Parent material for IGFEM * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "igfem_internal_field.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_IGFEM_HH_ #define AKANTU_MATERIAL_IGFEM_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelIGFEM; } namespace akantu { class MaterialIGFEM : public virtual Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEEngineTemplate MyFEEngineIGFEMType; MaterialIGFEM(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialIGFEM(); protected: /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void computeAllStresses(GhostType ghost_type = _not_ghost); virtual void extrapolateInternal(const ID & id, const Element & element, const Matrix & point, Matrix & extrapolated); /// apply a constant eigengrad_u everywhere in the material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & sub_mat_name, const GhostType = _not_ghost); /* ------------------------------------------------------------------------ */ /* MeshEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ virtual void computeQuadraturePointsCoordinates( ElementTypeMapArray & quadrature_points_coordinates, GhostType ghost_type) const; // virtual void onElementsAdded(const Array & element_list, // const NewElementsEvent & event) {}; // virtual void onElementsRemoved(const Array & element_list, - // const ElementTypeMapArray & + // const ElementTypeMapArray & // new_numbering, // const RemovedElementsEvent & event) {}; protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} void initialize(); template void setSubMaterial(const Array & element_list, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: const UInt nb_sub_materials; /// pointer to the solid mechanics model for igfem elements SolidMechanicsModelIGFEM * model; /// internal field of bool to know to which sub-material a quad point belongs IGFEMInternalField sub_material; /// material name of first sub-material std::string name_sub_mat_1; /// material name of first sub-material std::string name_sub_mat_2; /// map the index of the sub-materials to the names std::map sub_material_names; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_igfem_inline_impl.hh" } // namespace akantu #include "igfem_internal_field_tmpl.hh" #endif /* AKANTU_MATERIAL_IGFEM_HH_ */ diff --git a/extra_packages/igfem/src/material_igfem/material_igfem_elastic.cc b/extra_packages/igfem/src/material_igfem/material_igfem_elastic.cc index 085a2a94c..a331e82ab 100644 --- a/extra_packages/igfem/src/material_igfem/material_igfem_elastic.cc +++ b/extra_packages/igfem/src/material_igfem/material_igfem_elastic.cc @@ -1,241 +1,241 @@ /** * @file material_igfem_elastic.cc * * @author Aurelia Isabel Cuba Ramos * * * @brief Specializaton of material class for the igfem elastic material * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_igfem_elastic.hh" #include "material_elastic.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialIGFEMElastic::MaterialIGFEMElastic(SolidMechanicsModel & model, const ID & id) : Material(model, id), Parent(model, id), lambda("lambda", *this), mu("mu", *this), kpa("kappa", *this) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::initialize() { this->lambda.initialize(1); this->mu.initialize(1); this->kpa.initialize(1); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); /// insert the sub_material names into the map this->sub_material_names[0] = this->name_sub_mat_1; this->sub_material_names[1] = this->name_sub_mat_2; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::updateElasticInternals( const Array & element_list) { /// compute the Lamé constants for both sub-materials Vector lambda_per_sub_mat(this->nb_sub_materials); Vector mu_per_sub_mat(this->nb_sub_materials); Vector kpa_per_sub_mat(this->nb_sub_materials); for (UInt i = 0; i < this->nb_sub_materials; ++i) { ID mat_name = this->sub_material_names[i]; const MaterialElastic & mat = dynamic_cast &>( this->model->getMaterial(mat_name)); lambda_per_sub_mat(i) = mat.getLambda(); mu_per_sub_mat(i) = mat.getMu(); kpa_per_sub_mat(i) = mat.getKappa(); } for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType ghost_type = *g; /// loop over all types in the material - typedef ElementTypeMapArray::type_iterator iterator; + typedef ElementTypeMapArray::type_iterator iterator; iterator it = this->element_filter.firstType(spatial_dimension, ghost_type, _ek_igfem); iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type, _ek_igfem); /// loop over all types in the filter for (; it != last_type; ++it) { ElementType el_type = *it; if (el_type == _igfem_triangle_4) this->template setSubMaterial<_igfem_triangle_4>(element_list, ghost_type); else if (el_type == _igfem_triangle_5) this->template setSubMaterial<_igfem_triangle_5>(element_list, ghost_type); else AKANTU_ERROR("There is currently no other IGFEM type implemented"); UInt nb_element = this->element_filter(el_type, ghost_type).getSize(); UInt nb_quads = this->fem->getNbIntegrationPoints(el_type); /// get pointer to internals for given type Real * lambda_ptr = this->lambda(el_type, ghost_type).storage(); Real * mu_ptr = this->mu(el_type, ghost_type).storage(); Real * kpa_ptr = this->kpa(el_type, ghost_type).storage(); UInt * sub_mat_ptr = this->sub_material(el_type, ghost_type).storage(); for (UInt q = 0; q < nb_element * nb_quads; ++q, ++lambda_ptr, ++mu_ptr, ++kpa_ptr, ++sub_mat_ptr) { UInt index = *sub_mat_ptr; *lambda_ptr = lambda_per_sub_mat(index); *mu_ptr = mu_per_sub_mat(index); *kpa_ptr = kpa_per_sub_mat(index); } } } } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Parent::computeStress(el_type, ghost_type); if (!this->finite_deformation) { /// get pointer to internals Real * lambda_ptr = this->lambda(el_type, ghost_type).storage(); Real * mu_ptr = this->mu(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); this->computeStressOnQuad(grad_u, sigma, *lambda_ptr, *mu_ptr); ++lambda_ptr; ++mu_ptr; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } else { AKANTU_DEBUG_TO_IMPLEMENT(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::computeTangentModuli( __attribute__((unused)) ElementType el_type, Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type) { AKANTU_DEBUG_IN(); /// get pointer to internals Real * lambda_ptr = this->lambda(el_type, ghost_type).storage(); Real * mu_ptr = this->mu(el_type, ghost_type).storage(); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); this->computeTangentModuliOnQuad(tangent, *lambda_ptr, *mu_ptr); ++lambda_ptr; ++mu_ptr; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::computePotentialEnergy( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); // MaterialThermal::computePotentialEnergy(el_type, // ghost_type); // if(ghost_type != _not_ghost) return; // Array::scalar_iterator epot = this->potential_energy(el_type, // ghost_type).begin(); // if (!this->finite_deformation) { // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); // this->computePotentialEnergyOnQuad(grad_u, sigma, *epot); // ++epot; // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; // } else { // Matrix E(spatial_dimension, spatial_dimension); // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); // this->template gradUToGreenStrain(grad_u, E); // this->computePotentialEnergyOnQuad(E, sigma, *epot); // ++epot; // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; // } AKANTU_DEBUG_TO_IMPLEMENT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::computePotentialEnergyByElement( ElementType type, UInt index, Vector & epot_on_quad_points) { // Array::matrix_iterator gradu_it = // this->gradu(type).begin(spatial_dimension, // spatial_dimension); // Array::matrix_iterator gradu_end = // this->gradu(type).begin(spatial_dimension, // spatial_dimension); // Array::matrix_iterator stress_it = // this->stress(type).begin(spatial_dimension, // spatial_dimension); // if (this->finite_deformation) // stress_it = this->piola_kirchhoff_2(type).begin(spatial_dimension, // spatial_dimension); // UInt nb_quadrature_points = // this->model->getFEEngine().getNbQuadraturePoints(type); // gradu_it += index*nb_quadrature_points; // gradu_end += (index+1)*nb_quadrature_points; // stress_it += index*nb_quadrature_points; // Real * epot_quad = epot_on_quad_points.storage(); // Matrix grad_u(spatial_dimension, spatial_dimension); // for(;gradu_it != gradu_end; ++gradu_it, ++stress_it, ++epot_quad) { // if (this->finite_deformation) // this->template gradUToGreenStrain(*gradu_it, // grad_u); // else // grad_u.copy(*gradu_it); // this->computePotentialEnergyOnQuad(grad_u, *stress_it, *epot_quad); // } AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { Parent::onElementsAdded(element_list, event); updateElasticInternals(element_list); }; INSTANTIATE_MATERIAL(MaterialIGFEMElastic); } // namespace akantu diff --git a/extra_packages/igfem/src/material_igfem/material_igfem_iterative_stiffness_reduction.cc b/extra_packages/igfem/src/material_igfem/material_igfem_iterative_stiffness_reduction.cc index 34e564db6..7ab8bcd37 100644 --- a/extra_packages/igfem/src/material_igfem/material_igfem_iterative_stiffness_reduction.cc +++ b/extra_packages/igfem/src/material_igfem/material_igfem_iterative_stiffness_reduction.cc @@ -1,288 +1,288 @@ /** * @file material_igfem_iterative_stiffness_reduction.cc * @author Aurelia Isabel Cuba Ramos * @date Thu Mar 10 08:37:43 2016 * * @brief Implementation of igfem material iterative stiffness reduction * * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_igfem_iterative_stiffness_reduction.hh" #include namespace akantu { template /* -------------------------------------------------------------------------- */ MaterialIGFEMIterativeStiffnessReduction:: MaterialIGFEMIterativeStiffnessReduction(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialIGFEMSawToothDamage(model, id), eps_u("ultimate_strain", *this), reduction_step("damage_step", *this), D("tangent", *this), Gf(0.), crack_band_width(0.), max_reductions(0), reduction_constant(0.) { AKANTU_DEBUG_IN(); this->eps_u.initialize(1); this->D.initialize(1); this->reduction_step.initialize(1); this->internals_to_transfer.push_back("ultimate_strain"); this->internals_to_transfer.push_back("tangent"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMIterativeStiffnessReduction::initMaterial() { AKANTU_DEBUG_IN(); MaterialIGFEMSawToothDamage::initMaterial(); /// get the parameters for the sub-material that can be damaged ID mat_name = this->sub_material_names[1]; const Material & mat = this->model->getMaterial(mat_name); this->crack_band_width = mat.getParam("crack_band_width"); this->max_reductions = mat.getParam("max_reductions"); this->reduction_constant = mat.getParam("reduction_constant"); this->Gf = mat.getParam("Gf"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMIterativeStiffnessReduction:: computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// storage for the current stress Matrix sigma(spatial_dimension, spatial_dimension); /// Vector to store eigenvalues of current stress tensor Vector eigenvalues(spatial_dimension); /// iterators on the needed internal fields Array::const_scalar_iterator Sc_it = this->Sc(el_type, ghost_type).begin(); Array::scalar_iterator dam_it = this->damage(el_type, ghost_type).begin(); Array::scalar_iterator equivalent_stress_it = this->equivalent_stress(el_type, ghost_type).begin(); Array::const_matrix_iterator grad_u_it = grad_u.begin(spatial_dimension, spatial_dimension); Array::const_matrix_iterator grad_u_end = grad_u.end(spatial_dimension, spatial_dimension); Real * mu_ptr = this->mu(el_type, ghost_type).storage(); Real * lambda_ptr = this->lambda(el_type, ghost_type).storage(); /// loop over all the quadrature points and compute the equivalent stress for (; grad_u_it != grad_u_end; ++grad_u_it) { /// compute the stress sigma.zero(); MaterialIGFEMElastic::computeStressOnQuad( *grad_u_it, sigma, *lambda_ptr, *mu_ptr); MaterialIGFEMSawToothDamage< spatial_dimension>::computeDamageAndStressOnQuad(sigma, *dam_it); /// compute eigenvalues sigma.eig(eigenvalues); /// find max eigenvalue and normalize by tensile strength *equivalent_stress_it = *(std::max_element(eigenvalues.storage(), eigenvalues.storage() + spatial_dimension)) / (*Sc_it); ++Sc_it; ++equivalent_stress_it; ++dam_it; ++lambda_ptr; ++mu_ptr; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template UInt MaterialIGFEMIterativeStiffnessReduction< spatial_dimension>::updateDamage() { UInt nb_damaged_elements = 0; if (this->norm_max_equivalent_stress >= 1.) { AKANTU_DEBUG_IN(); /// update the damage only on non-ghosts elements! Doesn't make sense to /// update on ghost. GhostType ghost_type = _not_ghost; ; Mesh::type_iterator it = this->model->getFEEngine().getMesh().firstType( spatial_dimension, ghost_type, _ek_igfem); Mesh::type_iterator last_type = this->model->getFEEngine().getMesh().lastType(spatial_dimension, ghost_type, _ek_igfem); /// get the Young's modulus of the damageable sub-material ID mat_name = this->sub_material_names[1]; Real E = this->model->getMaterial(mat_name).template getParam("E"); /// loop over all the elements for (; it != last_type; ++it) { ElementType el_type = *it; /// get iterators on the needed internal fields - const Array & sub_mat = this->sub_material(el_type, ghost_type); - Array::const_scalar_iterator sub_mat_it = sub_mat.begin(); + const Array & sub_mat = this->sub_material(el_type, ghost_type); + Array::const_scalar_iterator sub_mat_it = sub_mat.begin(); Array::const_scalar_iterator equivalent_stress_it = this->equivalent_stress(el_type, ghost_type).begin(); Array::const_scalar_iterator equivalent_stress_end = this->equivalent_stress(el_type, ghost_type).end(); Array::scalar_iterator dam_it = this->damage(el_type, ghost_type).begin(); - Array::scalar_iterator reduction_it = + Array::scalar_iterator reduction_it = this->reduction_step(el_type, ghost_type).begin(); Array::const_scalar_iterator eps_u_it = this->eps_u(el_type, ghost_type).begin(); Array::scalar_iterator Sc_it = this->Sc(el_type, ghost_type).begin(); Array::const_scalar_iterator D_it = this->D(el_type, ghost_type).begin(); /// loop over all the elements of the given type UInt nb_element = this->element_filter(el_type, ghost_type).getSize(); UInt nb_quads = this->fem->getNbIntegrationPoints(el_type, ghost_type); bool damage_element = false; for (UInt e = 0; e < nb_element; ++e) { damage_element = false; /// check if damage occurs in the element for (UInt q = 0; q < nb_quads; ++q, ++reduction_it, ++sub_mat_it, ++equivalent_stress_it) { if (*equivalent_stress_it >= (1 - this->dam_tolerance) * this->norm_max_equivalent_stress && *sub_mat_it != 0) { /// check if this element can still be damaged if (*reduction_it == this->max_reductions) continue; damage_element = true; } } if (damage_element) { /// damage the element nb_damaged_elements += 1; sub_mat_it -= nb_quads; reduction_it -= nb_quads; for (UInt q = 0; q < nb_quads; ++q) { if (*sub_mat_it) { /// increment the counter of stiffness reduction steps *reduction_it += 1; if (*reduction_it == this->max_reductions) *dam_it = this->max_damage; else { /// update the damage on this quad *dam_it = 1. - (1. / std::pow(this->reduction_constant, *reduction_it)); /// update the stiffness on this quad *Sc_it = (*eps_u_it) * (1. - (*dam_it)) * E * (*D_it) / ((1. - (*dam_it)) * E + (*D_it)); } } ++sub_mat_it; ++dam_it; ++reduction_it; ++eps_u_it; ++Sc_it; ++D_it; } } else { dam_it += nb_quads; eps_u_it += nb_quads; Sc_it += nb_quads; D_it += nb_quads; } } } } StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&nb_damaged_elements, 1, _so_sum); AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMIterativeStiffnessReduction< spatial_dimension>::onElementsAdded(__attribute__((unused)) const Array & element_list, __attribute__((unused)) const NewElementsEvent & event) { MaterialIGFEMSawToothDamage::onElementsAdded(element_list, event); /// set the correct damage iteration step (is UInt->cannot be interpolated) Real val = 0.; for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType ghost_type = *g; /// loop over all types in the material - typedef ElementTypeMapArray::type_iterator iterator; + typedef ElementTypeMapArray::type_iterator iterator; iterator it = this->element_filter.firstType(spatial_dimension, ghost_type, _ek_igfem); iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type, _ek_igfem); /// loop over all types in the filter for (; it != last_type; ++it) { const ElementType el_type = *it; Array::scalar_iterator dam_it = this->damage(el_type, ghost_type).begin(); - Array::scalar_iterator reduction_it = + Array::scalar_iterator reduction_it = this->reduction_step(el_type, ghost_type).begin(); UInt nb_element = this->element_filter(el_type, ghost_type).getSize(); UInt nb_quads = this->fem->getNbIntegrationPoints(el_type); UInt * sub_mat_ptr = this->sub_material(el_type, ghost_type).storage(); for (UInt q = 0; q < nb_element * nb_quads; ++q, ++sub_mat_ptr, ++dam_it, ++reduction_it) { if (*sub_mat_ptr) { if (Math::are_float_equal(*dam_it, this->max_damage)) *reduction_it = this->max_reductions; else { for (UInt i = 0; i < this->max_reductions; ++i) { val = 1 - (1. / std::pow(this->reduction_constant, i)); if (Math::are_float_equal(val, *dam_it)) *reduction_it = i; } } } } } } } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialIGFEMIterativeStiffnessReduction); } // namespace akantu diff --git a/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage.cc b/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage.cc index 5d35bddf3..5936ca6b1 100644 --- a/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage.cc +++ b/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage.cc @@ -1,466 +1,466 @@ /** * @file material_igfem_saw_tooth_damage.cc * * @author Aurelia Isabel Cuba Ramos * * * @brief Implementation of the squentially linear saw-tooth damage model for * IGFEM elements * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_igfem_saw_tooth_damage.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialIGFEMSawToothDamage::MaterialIGFEMSawToothDamage( SolidMechanicsModel & model, const ID & id) : Material(model, id), Parent(model, id), Sc("Sc", *this), equivalent_stress("equivalent_stress", *this), norm_max_equivalent_stress(0) { AKANTU_DEBUG_IN(); this->Sc.initialize(1); this->equivalent_stress.initialize(1); this->damage.setElementKind(_ek_igfem); this->damage.setFEEngine(*(this->fem)); this->internals_to_transfer.push_back("damage"); this->internals_to_transfer.push_back("Sc"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMSawToothDamage::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); /// get the parameters for the sub-material that can be damaged ID mat_name = this->sub_material_names[1]; const Material & mat = this->model->getMaterial(mat_name); this->prescribed_dam = mat.getParam("prescribed_dam"); this->dam_threshold = mat.getParam("dam_threshold"); this->dam_tolerance = mat.getParam("dam_tolerance"); this->max_damage = mat.getParam("max_damage"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMSawToothDamage:: computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// Vector to store eigenvalues of current stress tensor Vector eigenvalues(spatial_dimension); Array::const_iterator Sc_it = Sc(el_type, ghost_type).begin(); Array::iterator equivalent_stress_it = equivalent_stress(el_type, ghost_type).begin(); Array::const_matrix_iterator grad_u_it = grad_u.begin(spatial_dimension, spatial_dimension); Array::const_matrix_iterator grad_u_end = grad_u.end(spatial_dimension, spatial_dimension); /// get pointer to internals Real * lambda_ptr = this->lambda(el_type, ghost_type).storage(); Real * mu_ptr = this->mu(el_type, ghost_type).storage(); Real * dam = this->damage(el_type, ghost_type).storage(); Matrix sigma(spatial_dimension, spatial_dimension); for (; grad_u_it != grad_u_end; ++grad_u_it) { MaterialIGFEMElastic::computeStressOnQuad( *grad_u_it, sigma, *lambda_ptr, *mu_ptr); computeDamageAndStressOnQuad(sigma, *dam); /// compute eigenvalues sigma.eig(eigenvalues); /// find max eigenvalue and normalize by tensile strength *equivalent_stress_it = *(std::max_element(eigenvalues.storage(), eigenvalues.storage() + spatial_dimension)) / *(Sc_it); ++Sc_it; ++equivalent_stress_it; ++dam; ++lambda_ptr; ++mu_ptr; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMSawToothDamage::computeAllStresses( GhostType ghost_type) { AKANTU_DEBUG_IN(); /// reset normalized maximum equivalent stress if (ghost_type == _not_ghost) norm_max_equivalent_stress = 0; Parent::computeAllStresses(ghost_type); /// find global Gauss point with highest stress StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&norm_max_equivalent_stress, 1, _so_max); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMSawToothDamage:: findMaxNormalizedEquivalentStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (ghost_type == _not_ghost) { const Array & e_stress = equivalent_stress(el_type); Array::const_iterator equivalent_stress_it = e_stress.begin(); Array::const_iterator equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); Array::iterator dam_it = dam.begin(); - Array & sub_mat = this->sub_material(el_type); - Array::iterator sub_mat_it = sub_mat.begin(); + Array & sub_mat = this->sub_material(el_type); + Array::iterator sub_mat_it = sub_mat.begin(); for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it, ++sub_mat_it) { /// check if max equivalent stress for this element type is greater than /// the current norm_max_eq_stress and if the element is not already fully /// damaged if ((*equivalent_stress_it > norm_max_equivalent_stress && *dam_it < max_damage) && *sub_mat_it != 0) { norm_max_equivalent_stress = *equivalent_stress_it; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMSawToothDamage::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Parent::computeStress(el_type, ghost_type); Real * dam = this->damage(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computeDamageAndStressOnQuad(sigma, *dam); ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; computeNormalizedEquivalentStress(this->gradu(el_type, ghost_type), el_type, ghost_type); norm_max_equivalent_stress = 0; findMaxNormalizedEquivalentStress(el_type, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template UInt MaterialIGFEMSawToothDamage::updateDamage() { UInt nb_damaged_elements = 0; AKANTU_DEBUG_ASSERT(prescribed_dam > 0., "Your prescribed damage must be greater than zero"); if (norm_max_equivalent_stress >= 1.) { AKANTU_DEBUG_IN(); GhostType ghost_type = _not_ghost; ; Mesh::type_iterator it = this->model->getMesh().firstType( spatial_dimension, ghost_type, _ek_igfem); Mesh::type_iterator last_type = this->model->getMesh().lastType( spatial_dimension, ghost_type, _ek_igfem); for (; it != last_type; ++it) { ElementType el_type = *it; - Array & sub_mat = this->sub_material(el_type); - Array::iterator sub_mat_it = sub_mat.begin(); + Array & sub_mat = this->sub_material(el_type); + Array::iterator sub_mat_it = sub_mat.begin(); const Array & e_stress = equivalent_stress(el_type); Array::const_iterator equivalent_stress_it = e_stress.begin(); Array::const_iterator equivalent_stress_end = e_stress.end(); Array & dam = this->damage(el_type); Array::iterator dam_it = dam.begin(); /// loop over all the elements of the given type UInt nb_element = this->element_filter(el_type, ghost_type).getSize(); UInt nb_quads = this->fem->getNbIntegrationPoints(el_type, ghost_type); bool damage_element = false; for (UInt e = 0; e < nb_element; ++e) { damage_element = false; /// check if damage occurs in the element for (UInt q = 0; q < nb_quads; ++q) { if (*equivalent_stress_it >= (1 - dam_tolerance) * norm_max_equivalent_stress && *sub_mat_it != 0) { damage_element = true; } ++sub_mat_it; ++equivalent_stress_it; } if (damage_element) { nb_damaged_elements += 1; /// damage the element sub_mat_it -= nb_quads; for (UInt q = 0; q < nb_quads; ++q) { if (*sub_mat_it) { if (*dam_it < dam_threshold) *dam_it += prescribed_dam; else *dam_it = max_damage; } ++sub_mat_it; ++dam_it; } } else { dam_it += nb_quads; } } } } StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); comm.allReduce(&nb_damaged_elements, 1, _so_sum); AKANTU_DEBUG_OUT(); return nb_damaged_elements; } /* -------------------------------------------------------------------------- */ // template // void // MaterialIGFEMSawToothDamage::updateEnergiesAfterDamage(ElementType // el_type, GhostType ghost_type) { // MaterialDamage::updateEnergies(el_type, ghost_type); // } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMSawToothDamage::onElementsAdded( __attribute__((unused)) const Array & element_list, __attribute__((unused)) const NewElementsEvent & event) { /// update elastic constants MaterialIGFEMElastic::onElementsAdded(element_list, event); IGFEMInternalField quadrature_points_coordinates( "quadrature_points_coordinates", *this); quadrature_points_coordinates.initialize(spatial_dimension); this->computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); this->computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); Array::const_iterator el_begin = element_list.begin(); Array::const_iterator el_end = element_list.end(); Element el; el.kind = _ek_igfem; /// loop over all the elements in the filter for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType ghost_type = *g; /// loop over all types in the material - typedef ElementTypeMapArray::type_iterator iterator; + typedef ElementTypeMapArray::type_iterator iterator; iterator it = this->element_filter.firstType(spatial_dimension, ghost_type, _ek_igfem); iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type, _ek_igfem); for (; it != last_type; ++it) { /// store the elements added to this material std::vector new_elements; - Array & elem_filter = this->element_filter(*it, ghost_type); + Array & elem_filter = this->element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quads = quadrature_points_coordinates.getFEEngine().getNbIntegrationPoints( *it, ghost_type); Array added_quads(0, spatial_dimension); const Array & quads = quadrature_points_coordinates(*it, ghost_type); Array::const_vector_iterator quad = quads.begin(spatial_dimension); el.type = *it; el.ghost_type = ghost_type; for (UInt e = 0; e < nb_element; ++e) { /// global number of element el.element = elem_filter(e); if (std::find(el_begin, el_end, el) != el_end) { new_elements.push_back(el); for (UInt q = 0; q < nb_quads; ++q) { added_quads.push_back(*quad); ++quad; } } else quad += nb_quads; } /// get the extrapolated values for (UInt i = 0; i < this->internals_to_transfer.size(); ++i) { if (!new_elements.size()) continue; const ID name = this->getID() + ":" + this->internals_to_transfer[i]; Array & internal = (*(this->internal_vectors_real[name]))(*it, ghost_type); UInt nb_component = internal.getNbComponent(); /// Array::vector_iterator internal_it = /// internal.begin(nb_component); Array extrapolated_internal_values(new_elements.size() * nb_quads, nb_component); this->model->transferInternalValues(this->internals_to_transfer[i], new_elements, added_quads, extrapolated_internal_values); UInt * sub_mat_ptr = this->sub_material(*it, ghost_type).storage(); Array::const_vector_iterator extrapolated_it = extrapolated_internal_values.begin(nb_component); Array::vector_iterator internal_it = internal.begin(nb_component); /// copy back the results for (UInt j = 0; j < new_elements.size(); ++j) { Element element_local = this->convertToLocalElement(new_elements[j]); for (UInt q = 0; q < nb_quads; ++q, ++extrapolated_it) { if (sub_mat_ptr[element_local.element * nb_quads + q]) internal_it[element_local.element * nb_quads + q] = *extrapolated_it; } } } } } } /* -------------------------------------------------------------------------- */ // template // void // MaterialIGFEMSawToothDamage::transferInternals(Material & // old_mat, // std::vector & element_pairs) // { // Element new_el_global; // Element new_el_local; // Element old_el_global; // Element old_el_local; // /// get the fe-engine of the old material // FEEngine & fem_old_mat = old_mat.getFEEngine(); // for (UInt e = 0; e < element_pairs.size(); ++e) { // new_el_global = element_pairs[e].first; // old_el_global = element_pairs[e].second; // /// get the number of the elements in their materials // old_el_local = old_el_global; -// Array & mat_local_numbering = +// Array & mat_local_numbering = // this->model->getMaterialLocalNumbering(old_el_global.type, // old_el_global.ghost_type); // old_el_local.element = mat_local_numbering(old_el_global.element); // new_el_local = this->convertToLocalElement(new_el_global); // UInt nb_old_quads = fem_old_mat.getNbIntegrationPoints(old_el_global.type, // old_el_global.ghost_type); // UInt nb_new_quads = this->fem->getNbIntegrationPoints(new_el_global.type, // new_el_global.ghost_type); // if (old_mat.isInternal("damage", Mesh::getKind(old_el_global.type)) // && old_mat.isInternal("Sc", Mesh::getKind(old_el_global.type)) ) { // UInt quad; // Vector el_old_damage(nb_old_quads); // Vector el_old_Sc(nb_old_quads); // Vector el_new_damage(nb_new_quads); // Vector el_new_Sc(nb_new_quads); // const Array & old_Sc = old_mat.getArray("Sc", old_el_global.type, // old_el_global.ghost_type); // const Array & old_damage = old_mat.getArray("damage", // old_el_global.type, old_el_global.ghost_type); // Array & new_Sc = this->Sc(new_el_global.type, // new_el_global.ghost_type); // Array & new_damage = this->damage(new_el_global.type, // new_el_global.ghost_type); // for (UInt q = 0; q < nb_old_quads; ++q) { // quad = old_el_local.element * nb_old_quads + q; // el_old_damage(q) = old_damage(quad); // el_old_Sc(q) = old_Sc(quad); // } // this->interpolateInternal(new_el_global, old_el_global, el_new_damage, // el_old_damage, nb_new_quads, nb_old_quads); // this->interpolateInternal(new_el_global, old_el_global, el_new_Sc, // el_old_Sc, nb_new_quads, nb_old_quads); // for (UInt q = 0; q < nb_new_quads; ++q) { // quad = new_el_local.element * nb_new_quads + q; // if (this->sub_material(new_el_global.type,new_el_global.ghost_type)(quad)) { // new_damage(quad) = el_new_damage(q); // new_Sc(quad) = el_new_damage(q); // } // } // } // else // AKANTU_DEBUG_ASSERT((!old_mat.isInternal("damage", // Mesh::getKind(old_el_global.type)) // && !old_mat.isInternal("Sc", Mesh::getKind(old_el_global.type)) // ), // "old material has damage or Sc but not both!!!!"); // } // } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialIGFEMSawToothDamage); } // namespace akantu /* /// get the material index in the model from this current material UInt this_mat_idx = this->model->getMaterialIndex(this->name); /// number of elements in this event UInt nb_new_elements = element_list.getSize(); // const Array & old_elements = igfem_event->getOldElementsList(); // std::map > elements_by_old_mat; /// store the old elements sorted by their material for (UInt e = 0; e < nb_new_elements; ++e) { const Element new_el = element_list(e); - const Array & mat_idx = + const Array & mat_idx = this->model->getMaterialByElement(new_el.type, new_el.ghost_type); if ( mat_idx(new_el.element) != this_mat_idx ) continue; /// new element is not part of this material: nothing to be done /// get the corresponding old element and store new and old one as pair const Element old_el = old_elements(e); ElementPair el_pair(new_el, old_el); - const Array & old_mat_idx = + const Array & old_mat_idx = this->model->getMaterialByElement(old_el.type, old_el.ghost_type); UInt mat_old_idx = old_mat_idx(old_el.element); elements_by_old_mat[mat_old_idx].push_back(el_pair); } /// loop over all the element pairs in the map for (std::map >::iterator map_it = elements_by_old_mat.begin(); map_it != elements_by_old_mat.end(); ++map_it) { /// get the vector of old and new element pairs std::vector & element_pairs = map_it->second; Material & old_mat = this->model->getMaterial(map_it->first); this->transferInternals(old_mat, element_pairs); } */ diff --git a/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage_inline_impl.hh b/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage_inline_impl.hh index 45d8c6bac..8bc139cb9 100644 --- a/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage_inline_impl.hh +++ b/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage_inline_impl.hh @@ -1,66 +1,66 @@ /** * @file material_igfem_saw_tooth_damage_inline_impl.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Implementation of inline functions of the squentially linear * saw-tooth damage model for IGFEM elements * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ template inline void MaterialIGFEMSawToothDamage::computeDamageAndStressOnQuad( Matrix & sigma, Real & dam) { sigma *= 1 - dam; } /* -------------------------------------------------------------------------- */ template UInt MaterialIGFEMSawToothDamage::updateDamage( UInt quad_index, const Real eq_stress, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_ASSERT(prescribed_dam > 0., "Your prescribed damage must be greater than zero"); Array & dam = this->damage(el_type, ghost_type); Real & dam_on_quad = dam(quad_index); /// check if damage occurs if (equivalent_stress(el_type, ghost_type)(quad_index) >= (1 - dam_tolerance) * norm_max_equivalent_stress) { /// damage the entire sub-element -> get the element index UInt el_index = quad_index / this->element_filter(el_type, ghost_type).getSize(); UInt nb_quads = this->fem->getNbIntegrationPoints(el_type, ghost_type); UInt start_idx = el_index * nb_quads; - Array & sub_mat = this->sub_material(el_type, ghost_type); + Array & sub_mat = this->sub_material(el_type, ghost_type); UInt damaged_quads = 0; if (dam_on_quad < dam_threshold) { for (UInt q = 0; q < nb_quads; ++q, ++start_idx) { if (sub_mat(start_idx)) { dam(start_idx) += prescribed_dam; damaged_quads += 1; } } } else { for (UInt q = 0; q < nb_quads; ++q, ++start_idx) { if (sub_mat(start_idx)) { dam(start_idx) += max_damage; damaged_quads += 1; } } } return damaged_quads; } return 0; } /* -------------------------------------------------------------------------- */ diff --git a/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel.hh b/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel.hh index 10a101e6e..c052ed89e 100644 --- a/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel.hh +++ b/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel.hh @@ -1,220 +1,220 @@ /** * @file mesh_igfem_spherical_growing_gel.hh * * @author Clement Roux-Langlois * * @date creation: Mon Jul 13 2015 * * @brief Computation of mesh intersection with sphere(s) and growing of these * spheres. This class handle the intersectors templated for every * element * types. * * * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ //#if 0 #ifndef AKANTU_MESH_IGFEM_SPHERICAL_GROWING_GEL_HH_ #define AKANTU_MESH_IGFEM_SPHERICAL_GROWING_GEL_HH_ #include "aka_common.hh" #include "mesh_sphere_intersector.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* classes for new igfem elements mesh events */ /* -------------------------------------------------------------------------- */ class NewIGFEMElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(OldElementsList, old_elements, Array &); AKANTU_GET_MACRO(OldElementsList, old_elements, const Array &); protected: Array old_elements; }; class NewIGFEMNodesEvent : public NewNodesEvent { public: - void setNewNodePerElem(const Array & new_node_per_elem) { + void setNewNodePerElem(const Array & new_node_per_elem) { this->new_node_per_elem = &new_node_per_elem; } void setType(ElementType new_type) { type = new_type; } void setGhostType(GhostType new_type) { ghost_type = new_type; } - AKANTU_GET_MACRO(NewNodePerElem, *new_node_per_elem, const Array &); + AKANTU_GET_MACRO(NewNodePerElem, *new_node_per_elem, const Array &); AKANTU_GET_MACRO(ElementType, type, ElementType); AKANTU_GET_MACRO(GhostType, ghost_type, GhostType); protected: ElementType type; GhostType ghost_type; - const Array * new_node_per_elem; + const Array * new_node_per_elem; }; /* -------------------------------------------------------------------------- */ template class MeshIgfemSphericalGrowingGel { // definition of the element list #define ELEMENT_LIST (_triangle_3)(_igfem_triangle_4)(_igfem_triangle_5) // Solution 1 /* #define INTERSECTOR_DEFINITION(type) \ MeshSphereIntersector<2, type> intersector##type(mesh);*/ // Solution 2 #define INSTANTIATOR(_type) \ intersectors(_type, ghost_type) = \ new MeshSphereIntersector(this->mesh) /* -------------------------------------------------------------------------- */ /* Constructor/Destructor */ /* -------------------------------------------------------------------------- */ public: /// Construct from mesh MeshIgfemSphericalGrowingGel(Mesh & mesh); /// Destructor ~MeshIgfemSphericalGrowingGel() { for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator it = mesh.firstType(dim, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(dim, ghost_type, _ek_not_defined); for (; it != end; ++it) { delete intersectors(*it, ghost_type); } } } /* -------------------------------------------------------------------------- */ /* Methods */ /* -------------------------------------------------------------------------- */ public: void init(); /// Remove the additionnal nodes void removeAdditionalNodes(); /// Compute the intersection points between the mesh and the query list for /// all element types and send the NewNodeEvent void computeMeshQueryListIntersectionPoint( const std::list & query_list); /// increase sphere radius and build the new intersection points between the /// mesh and the query list for all element types and send the NewNodeEvent void computeMeshQueryListIntersectionPoint( const std::list & query_list, Real inf_fact) { std::list::const_iterator query_it = query_list.begin(), query_end = query_list.end(); std::list sphere_list; for (; query_it != query_end; ++query_it) { SK::Sphere_3 sphere(query_it->center(), query_it->squared_radius() * inf_fact * inf_fact); sphere_list.push_back(sphere); } computeMeshQueryListIntersectionPoint(sphere_list); } /// Build the IGFEM mesh from intersection points void buildIGFEMMesh(); /// Build the IGFEM mesh from spheres void buildIGFEMMeshFromSpheres(const std::list & query_list) { computeMeshQueryListIntersectionPoint(query_list); buildIGFEMMesh(); } /// Build the IGFEM mesh from spheres with increase factor void buildIGFEMMeshFromSpheres(const std::list & query_list, Real inf_fact) { computeMeshQueryListIntersectionPoint(query_list, inf_fact); buildIGFEMMesh(); } /// set the distributed synchronizer void setDistributedSynchronizer(DistributedSynchronizer * dist) { synchronizer = dist; buildSegmentConnectivityToNodeType(); } /// update node type - void updateNodeType(const Array & nodes_list, - const Array & new_node_per_elem, ElementType type, + void updateNodeType(const Array & nodes_list, + const Array & new_node_per_elem, ElementType type, GhostType ghost_type); protected: /// Build the unordered_map needed to assign the node type to new nodes in /// parallel void buildSegmentConnectivityToNodeType(); /* -------------------------------------------------------------------------- */ /* Accessors */ /* -------------------------------------------------------------------------- */ public: AKANTU_GET_MACRO(NbStandardNodes, nb_nodes_fem, UInt); AKANTU_GET_MACRO(NbEnrichedNodes, nb_enriched_nodes, UInt); /* -------------------------------------------------------------------------- */ /* Class Members */ /* -------------------------------------------------------------------------- */ protected: /// Mesh used to construct the primitives Mesh & mesh; /// number of fem nodes in the initial mesh UInt nb_nodes_fem; /// number of enriched nodes before intersection UInt nb_enriched_nodes; // Solution 2 /// map of the elements types in the mesh and the corresponding intersectors ElementTypeMap *> intersectors; /// Map linking pairs of nodes to a node type. The pairs of nodes /// contain the connectivity of the primitive segments that are /// intersected. unordered_map, Int>::type segment_conn_to_node_type; /// Pointer to the distributed synchronizer of the model DistributedSynchronizer * synchronizer; }; } // namespace akantu #include "mesh_igfem_spherical_growing_gel_tmpl.hh" #endif // AKANTU_MESH_IGFEM_SPHERICAL_GROWING_GEL_HH_ //#endif // diff --git a/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh b/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh index d9e9ce065..09b591bfb 100644 --- a/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh +++ b/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh @@ -1,531 +1,531 @@ /** * @file mesh_igfem_spherical_growing_gel.cc * @author Marco Vocialta * @date Mon Sep 21 12:42:11 2015 * * @brief Implementation of the functions of MeshIgfemSphericalGrowingGel * * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_utils.hh" #include namespace akantu { /* -------------------------------------------------------------------------- */ template MeshIgfemSphericalGrowingGel::MeshIgfemSphericalGrowingGel(Mesh & mesh) : mesh(mesh), nb_nodes_fem(mesh.getNbNodes()), nb_enriched_nodes(0), synchronizer(NULL) {} /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::init() { nb_nodes_fem = mesh.getNbNodes(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator tit = mesh.firstType(dim, ghost_type); Mesh::type_iterator tend = mesh.lastType(dim, ghost_type); for (; tit != tend; ++tit) { // loop to add corresponding IGFEM element types if (*tit == _triangle_3) { mesh.addConnectivityType(_igfem_triangle_4, ghost_type); mesh.addConnectivityType(_igfem_triangle_5, ghost_type); } else AKANTU_ERROR("Not ready for mesh type " << *tit); } tit = mesh.firstType(dim, ghost_type, _ek_not_defined); tend = mesh.lastType(dim, ghost_type, _ek_not_defined); for (; tit != tend; ++tit) { AKANTU_BOOST_LIST_SWITCH(INSTANTIATOR, ELEMENT_LIST, *tit); } } } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::computeMeshQueryListIntersectionPoint( const std::list & query_list) { /// store number of currently enriched nodes this->nb_enriched_nodes = mesh.getNbNodes() - nb_nodes_fem; UInt nb_old_nodes = mesh.getNbNodes(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator iit = mesh.firstType(dim, ghost_type, _ek_not_defined); Mesh::type_iterator iend = mesh.lastType(dim, ghost_type, _ek_not_defined); for (; iit != iend; ++iit) { MeshAbstractIntersector & intersector = *intersectors(*iit, ghost_type); intersector.constructData(ghost_type); intersector.computeMeshQueryListIntersectionPoint(query_list, nb_old_nodes); const Array intersection_points_current_type = *(intersector.getIntersectionPoints()); - const Array & new_node_per_elem = intersector.getNewNodePerElem(); + const Array & new_node_per_elem = intersector.getNewNodePerElem(); /// Send the new node event UInt new_points = intersection_points_current_type.getSize(); StaticCommunicator::getStaticCommunicator().allReduce(&new_points, 1, _so_sum); if (new_points > 0) { Array & nodes = this->mesh.getNodes(); UInt nb_node = nodes.getSize(); Array::const_vector_iterator intersec_p_it = intersection_points_current_type.begin(dim); NewIGFEMNodesEvent new_nodes_event; for (UInt in = 0; in < intersection_points_current_type.getSize(); ++in, ++intersec_p_it) { nodes.push_back(*intersec_p_it); new_nodes_event.getList().push_back(nb_node + in); } new_nodes_event.setNewNodePerElem(new_node_per_elem); new_nodes_event.setType(*iit); new_nodes_event.setGhostType(ghost_type); this->mesh.sendEvent(new_nodes_event); } } } } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::removeAdditionalNodes() { AKANTU_DEBUG_IN(); UInt total_nodes = this->mesh.getNbNodes(); UInt nb_new_enriched_nodes = total_nodes - this->nb_enriched_nodes - this->nb_nodes_fem; UInt old_total_nodes = this->nb_nodes_fem + this->nb_enriched_nodes; UInt total_new_nodes = nb_new_enriched_nodes; StaticCommunicator::getStaticCommunicator().allReduce(&total_new_nodes, 1, _so_sum); if (total_new_nodes == 0) return; RemovedNodesEvent remove_nodes(this->mesh); - Array & nodes_removed = remove_nodes.getList(); - Array & new_numbering = remove_nodes.getNewNumbering(); + Array & nodes_removed = remove_nodes.getList(); + Array & new_numbering = remove_nodes.getNewNumbering(); for (UInt nnod = 0; nnod < this->nb_nodes_fem; ++nnod) { new_numbering(nnod) = nnod; } for (UInt nnod = nb_nodes_fem; nnod < old_total_nodes; ++nnod) { new_numbering(nnod) = UInt(-1); nodes_removed.push_back(nnod); } for (UInt nnod = 0; nnod < nb_new_enriched_nodes; ++nnod) { new_numbering(nnod + old_total_nodes) = this->nb_nodes_fem + nnod; } if (nb_new_enriched_nodes > 0) { for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType)gt; Mesh::type_iterator it = mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; - Array & connectivity_array = + Array & connectivity_array = mesh.getConnectivity(type, ghost_type); UInt nb_nodes_per_element = connectivity_array.getNbComponent(); - Array::vector_iterator conn_it = + Array::vector_iterator conn_it = connectivity_array.begin(nb_nodes_per_element); - Array::vector_iterator conn_end = + Array::vector_iterator conn_end = connectivity_array.end(nb_nodes_per_element); for (; conn_it != conn_end; ++conn_it) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { (*conn_it)(n) = new_numbering((*conn_it)(n)); } } } } } this->mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::buildIGFEMMesh() { AKANTU_DEBUG_IN(); NewIGFEMElementsEvent new_elements_event; UInt total_new_elements = 0; Array & new_elements_list = new_elements_event.getList(); Array & old_elements_list = new_elements_event.getOldElementsList(); RemovedElementsEvent removed_elements_event(this->mesh); ChangedElementsEvent changed_elements_event(this->mesh); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; mesh.initElementTypeMapArray(removed_elements_event.getNewNumbering(), 1, dim, ghost_type, false, _ek_not_defined); mesh.initElementTypeMapArray(changed_elements_event.getNewNumbering(), 1, dim, ghost_type, false, _ek_not_defined); /// loop over all types in the mesh for a given ghost type Mesh::type_iterator iit = mesh.firstType(dim, ghost_type, _ek_not_defined); Mesh::type_iterator iend = mesh.lastType(dim, ghost_type, _ek_not_defined); for (; iit != iend; ++iit) { ElementType type = *iit; MeshAbstractIntersector & intersector = *intersectors(type, ghost_type); - const Array & new_node_per_elem = intersector.getNewNodePerElem(); + const Array & new_node_per_elem = intersector.getNewNodePerElem(); UInt n_new_el = 0; - Array & connectivity = this->mesh.getConnectivity(type, ghost_type); + Array & connectivity = this->mesh.getConnectivity(type, ghost_type); /// get the connectivities of all types that that may transform - Array & connec_igfem_tri_4 = + Array & connec_igfem_tri_4 = this->mesh.getConnectivity(_igfem_triangle_4, ghost_type); - Array & connec_igfem_tri_5 = + Array & connec_igfem_tri_5 = this->mesh.getConnectivity(_igfem_triangle_5, ghost_type); - Array & connec_tri_3 = + Array & connec_tri_3 = this->mesh.getConnectivity(_triangle_3, ghost_type); /// create elements to store the newly generated elements Element el_tri_3(_triangle_3, 0, ghost_type, _ek_regular); Element el_igfem_tri_4(_igfem_triangle_4, 0, ghost_type, _ek_igfem); Element el_igfem_tri5(_igfem_triangle_5, 0, ghost_type, _ek_igfem); - Array & new_numbering = + Array & new_numbering = removed_elements_event.getNewNumbering(type, ghost_type); new_numbering.resize(connectivity.getSize()); /// container for element to be removed Element old_element(type, 0, ghost_type, Mesh::getKind(type)); for (UInt nel = 0; nel != new_node_per_elem.getSize(); ++nel) { /// a former IGFEM triangle is transformed into a regular triangle if ((type != _triangle_3) && (new_node_per_elem(nel, 0) == 0)) { Vector connec_new_elem(3); connec_new_elem(0) = connectivity(nel, 0); connec_new_elem(1) = connectivity(nel, 1); connec_new_elem(2) = connectivity(nel, 2); /// add the new element in the mesh UInt nb_triangles_3 = connec_tri_3.getSize(); connec_tri_3.push_back(connec_new_elem); el_tri_3.element = nb_triangles_3; new_elements_list.push_back(el_tri_3); /// the number of the old element in the mesh old_element.element = nel; old_elements_list.push_back(old_element); new_numbering(nel) = UInt(-1); ++n_new_el; } /// element is enriched else if (new_node_per_elem(nel, 0) != 0) { switch (new_node_per_elem(nel, 0)) { /// new element is of type igfem_triangle_4 case 1: { Vector connec_new_elem(4); switch (new_node_per_elem(nel, 2)) { case 0: connec_new_elem(0) = connectivity(nel, 2); connec_new_elem(1) = connectivity(nel, 0); connec_new_elem(2) = connectivity(nel, 1); break; case 1: connec_new_elem(0) = connectivity(nel, 0); connec_new_elem(1) = connectivity(nel, 1); connec_new_elem(2) = connectivity(nel, 2); break; case 2: connec_new_elem(0) = connectivity(nel, 1); connec_new_elem(1) = connectivity(nel, 2); connec_new_elem(2) = connectivity(nel, 0); break; default: AKANTU_ERROR("A triangle has only 3 segments not " << new_node_per_elem(nel, 2)); break; } connec_new_elem(3) = new_node_per_elem(nel, 1); UInt nb_igfem_triangles_4 = connec_igfem_tri_4.getSize(); connec_igfem_tri_4.push_back(connec_new_elem); el_igfem_tri_4.element = nb_igfem_triangles_4; new_elements_list.push_back(el_igfem_tri_4); break; } /// new element is of type igfem_triangle_5 case 2: { Vector connec_new_elem(5); if ((new_node_per_elem(nel, 2) == 0) && (new_node_per_elem(nel, 4) == 1)) { connec_new_elem(0) = connectivity(nel, 1); connec_new_elem(1) = connectivity(nel, 2); connec_new_elem(2) = connectivity(nel, 0); connec_new_elem(3) = new_node_per_elem(nel, 3); connec_new_elem(4) = new_node_per_elem(nel, 1); } else if ((new_node_per_elem(nel, 2) == 0) && (new_node_per_elem(nel, 4) == 2)) { connec_new_elem(0) = connectivity(nel, 0); connec_new_elem(1) = connectivity(nel, 1); connec_new_elem(2) = connectivity(nel, 2); connec_new_elem(3) = new_node_per_elem(nel, 1); connec_new_elem(4) = new_node_per_elem(nel, 3); } else if ((new_node_per_elem(nel, 2) == 1) && (new_node_per_elem(nel, 4) == 2)) { connec_new_elem(0) = connectivity(nel, 2); connec_new_elem(1) = connectivity(nel, 0); connec_new_elem(2) = connectivity(nel, 1); connec_new_elem(3) = new_node_per_elem(nel, 3); connec_new_elem(4) = new_node_per_elem(nel, 1); } else if ((new_node_per_elem(nel, 2) == 1) && (new_node_per_elem(nel, 4) == 0)) { connec_new_elem(0) = connectivity(nel, 1); connec_new_elem(1) = connectivity(nel, 2); connec_new_elem(2) = connectivity(nel, 0); connec_new_elem(3) = new_node_per_elem(nel, 1); connec_new_elem(4) = new_node_per_elem(nel, 3); } else if ((new_node_per_elem(nel, 2) == 2) && (new_node_per_elem(nel, 4) == 0)) { connec_new_elem(0) = connectivity(nel, 0); connec_new_elem(1) = connectivity(nel, 1); connec_new_elem(2) = connectivity(nel, 2); connec_new_elem(3) = new_node_per_elem(nel, 3); connec_new_elem(4) = new_node_per_elem(nel, 1); } else if ((new_node_per_elem(nel, 2) == 2) && (new_node_per_elem(nel, 4) == 1)) { connec_new_elem(0) = connectivity(nel, 2); connec_new_elem(1) = connectivity(nel, 0); connec_new_elem(2) = connectivity(nel, 1); connec_new_elem(3) = new_node_per_elem(nel, 1); connec_new_elem(4) = new_node_per_elem(nel, 3); } else AKANTU_ERROR("A triangle has only 3 segments (0 to 2) not " << new_node_per_elem(nel, 2) << "and" << new_node_per_elem(nel, 4)); UInt nb_igfem_triangles_5 = connec_igfem_tri_5.getSize(); connec_igfem_tri_5.push_back(connec_new_elem); el_igfem_tri5.element = nb_igfem_triangles_5; new_elements_list.push_back(el_igfem_tri5); break; } default: AKANTU_ERROR("IGFEM cannot add " << new_node_per_elem(nel, 0) << " nodes to triangles"); break; } old_element.element = nel; old_elements_list.push_back(old_element); new_numbering(nel) = UInt(-1); ++n_new_el; } } total_new_elements += n_new_el; } } /// update new numbering for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; /// loop over all types in the mesh for a given ghost type Mesh::type_iterator iit = mesh.firstType(dim, ghost_type, _ek_not_defined); Mesh::type_iterator iend = mesh.lastType(dim, ghost_type, _ek_not_defined); for (; iit != iend; ++iit) { ElementType type = *iit; - Array & new_numbering = + Array & new_numbering = removed_elements_event.getNewNumbering(type, ghost_type); UInt el_index = 0; UInt nb_element = this->mesh.getNbElement(type, ghost_type); new_numbering.resize(nb_element); for (UInt e = 0; e < nb_element; ++e) { if (new_numbering(e) != UInt(-1)) { new_numbering(e) = el_index; ++el_index; } } changed_elements_event.getNewNumbering(type, ghost_type) .copy(new_numbering); } } StaticCommunicator::getStaticCommunicator().allReduce(&total_new_elements, 1, _so_sum); if (total_new_elements > 0) { changed_elements_event.getListOld().copy( new_elements_event.getOldElementsList()); changed_elements_event.getListNew().copy(new_elements_event.getList()); this->mesh.sendEvent(changed_elements_event); this->mesh.sendEvent(new_elements_event); Array & removed_list = removed_elements_event.getList(); removed_list.copy(new_elements_event.getOldElementsList()); this->mesh.sendEvent(removed_elements_event); } removeAdditionalNodes(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::buildSegmentConnectivityToNodeType() { Mesh mesh_facets(mesh.initMeshFacets()); MeshUtils::buildSegmentToNodeType(mesh, mesh_facets, synchronizer); // only the ghost elements are considered for (UInt g = _not_ghost; g <= _ghost; ++g) { GhostType ghost_type = (GhostType)g; Mesh::type_iterator it = mesh_facets.firstType(1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(1, ghost_type); for (; it != end; ++it) { ElementType type = *it; const Array & segment_to_nodetype = mesh_facets.getData("segment_to_nodetype", type, ghost_type); - const Array & segment_connectivity = + const Array & segment_connectivity = mesh_facets.getConnectivity(type, ghost_type); // looping over all the segments Array::const_vector_iterator conn_it = segment_connectivity.begin(segment_connectivity.getNbComponent()); Array::const_vector_iterator conn_end = segment_connectivity.end(segment_connectivity.getNbComponent()); UInt seg_index = 0; UInt n[2]; for (; conn_it != conn_end; ++conn_it, ++seg_index) { Int seg_type = segment_to_nodetype(seg_index); n[0] = (*conn_it)(0); n[1] = (*conn_it)(1); if ((mesh.isMasterNode(n[0]) || mesh.isSlaveNode(n[0])) && (mesh.isMasterNode(n[1]) || mesh.isSlaveNode(n[1]))) { std::sort(n, n + 2); segment_conn_to_node_type[std::make_pair(n[0], n[1])] = seg_type; } } } } } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::updateNodeType( const Array & nodes_list, const Array & new_node_per_elem, ElementType type, GhostType ghost_type) { Array & nodes_type = mesh.getNodesType(); UInt old_nodes = nodes_type.getSize(); UInt new_nodes = nodes_list.getSize(); // exit this function if the simulation in run in serial if (old_nodes == 0 || new_nodes == 0) return; nodes_type.resize(old_nodes + new_nodes); Array checked_node(new_nodes, 1, false); UInt nb_prim_by_el = 0; if ((type == _triangle_3) || (type == _igfem_triangle_4) || (type == _igfem_triangle_5)) { nb_prim_by_el = 3; } else { AKANTU_ERROR("Not ready for mesh type " << type); } // determine the node type for the local, master and slave nodes const Array & connectivity = mesh.getConnectivity(type, ghost_type); unordered_map, Int>::type::iterator seg_type_it; unordered_map, Int>::type::iterator seg_type_end = segment_conn_to_node_type.end(); for (UInt el = 0; el < new_node_per_elem.getSize(); ++el) { UInt nb_nodes = new_node_per_elem(el, 0); for (UInt n = 0; n < nb_nodes; ++n) { UInt node_index = new_node_per_elem(el, 2 * n + 1); if (node_index < old_nodes || checked_node(node_index - old_nodes)) continue; // get the elemental connectivity of the segment associated to the node UInt segment_index = new_node_per_elem(el, 2 * n + 2); UInt extreme_nodes[2]; extreme_nodes[0] = segment_index; extreme_nodes[1] = segment_index + 1; if (extreme_nodes[1] == nb_prim_by_el) extreme_nodes[1] = 0; // get the connectivity of the segment extreme_nodes[0] = connectivity(el, extreme_nodes[0]); extreme_nodes[1] = connectivity(el, extreme_nodes[1]); // if one extreme nodes is pure ghost, then also the new node is pure // ghost if (mesh.isPureGhostNode(extreme_nodes[0]) || mesh.isPureGhostNode(extreme_nodes[1])) nodes_type(node_index) = -3; // if on of the two extreme nodes is local, then also the new node is // local else if (mesh.isLocalNode(extreme_nodes[0]) || mesh.isLocalNode(extreme_nodes[1])) nodes_type(node_index) = -1; // otherwise use the value stored in the map else { std::sort(extreme_nodes, extreme_nodes + 2); seg_type_it = segment_conn_to_node_type.find( std::make_pair(extreme_nodes[0], extreme_nodes[1])); AKANTU_DEBUG_ASSERT(seg_type_it != seg_type_end, "Segment not found"); nodes_type(node_index) = seg_type_it->second; } checked_node(node_index - old_nodes) = true; } } AKANTU_DEBUG_ASSERT(std::accumulate(checked_node.begin(), checked_node.end(), 0) == Int(checked_node.getSize()), "Not all new nodes were assigned a node type"); } } // namespace akantu diff --git a/packages/eigen.cmake b/packages/eigen.cmake index 841c85a63..96d65c88a 100644 --- a/packages/eigen.cmake +++ b/packages/eigen.cmake @@ -1,8 +1,7 @@ package_declare(Eigen3 EXTERNAL DESCRIPTION "Add Eigen3 dependency to akantu" EXTRA_PACKAGE_OPTIONS ARGS 3.3 COMPILE_FLAGS CXX -DEIGEN_MAX_CPP_VER=14 - COMPILE_FLAGS CXX -DEIGEN_DEFAULT_DENSE_INDEX_TYPE=akantu::Idx ) mark_as_advanced(Eigen3_DIR) diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh index a0f2de1a9..ef7835502 100644 --- a/src/common/aka_array_tmpl.hh +++ b/src/common/aka_array_tmpl.hh @@ -1,1041 +1,1044 @@ /** * @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_ /* -------------------------------------------------------------------------- */ #include "aka_array.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(UInt new_size, UInt 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.data()[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 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.data(), new_elem.size(), + std::copy_n(new_elem.derived().data(), new_elem.size(), values + this->nb_component * (this->size_ - 1)); } /// 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 UInt 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(Int 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(); } /* -------------------------------------------------------------------------- */ /** * 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 & other) { AKANTU_DEBUG_ASSERT((this->size_ == other.size_) && (this->nb_component == other.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 & other) { AKANTU_DEBUG_ASSERT((this->size_ == other.size()) && (this->nb_component == other.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 Array & Array::operator*=(const T & alpha) { 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 == array.data()) { return true; } return std::equal(this->values, this->values + this->size_ * this->nb_component, array.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