diff --git a/examples/contact_mechanics/CMakeLists.txt b/examples/contact_mechanics/CMakeLists.txt index be4a64e24..52c514b67 100644 --- a/examples/contact_mechanics/CMakeLists.txt +++ b/examples/contact_mechanics/CMakeLists.txt @@ -1,60 +1,60 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Mon Jan 18 2016 # # @brief configuration for heat transfer example # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_mesh(hertz_2d hertz_2d.geo 2 1) add_mesh(hertz_implicit_2d hertz_implicit_2d.geo 2 1) add_mesh(hertz_3d hertz_3d.geo 3 1) register_example(implicit_penalty_2d SOURCES implicit_penalty_2d.cc DEPENDS hertz_implicit_2d FILES_TO_COPY material_implicit.dat ) - register_example(explicit_penalty_2d SOURCES explicit_penalty_2d.cc DEPENDS hertz_2d FILES_TO_COPY material.dat ) register_example(explicit_penalty_3d SOURCES explicit_penalty_3d.cc DEPENDS hertz_3d FILES_TO_COPY material.dat ) register_example(solid_contact_explicit SOURCES solid_contact_explicit.cc DEPENDS hertz_2d FILES_TO_COPY material.dat ) + diff --git a/examples/contact_mechanics/solid_contact_explicit.cc b/examples/contact_mechanics/solid_contact_explicit.cc index 34b6ad944..b7e124d65 100644 --- a/examples/contact_mechanics/solid_contact_explicit.cc +++ b/examples/contact_mechanics/solid_contact_explicit.cc @@ -1,131 +1,200 @@ /** * @file contact_mechanics_penalty.cc * * @author Mohit Pundir * * @date creation: Mon Jan 24 2019 * @date last modification: Mon Jan 24 2019 * * @brief contact mechanics model with penalty resolution * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "non_linear_solver.hh" #include "contact_mechanics_model.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" +#include "solver_callback.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; + +class ContactSolverCallback + : public SolverCallback { + +public: + + ContactSolverCallback(SolidMechanicsModel &, ContactMechanicsModel &); + +public: + + void assembleMatrix(const ID &) override; + + void assembleResidual() override; + + void assembleLumpedMatrix(const ID &) override; + + MatrixType getMatrixType(const ID &) override; + +private: + + SolidMechanicsModel & solid; + + ContactMechanicsModel & contact; + +}; + +/* -------------------------------------------------------------------------- */ +ContactSolverCallback::ContactSolverCallback(SolidMechanicsModel & solid, + ContactMechanicsModel & contact) + : SolverCallback(), solid(solid), contact(contact) { + + +} + +/* -------------------------------------------------------------------------- */ +void ContactSolverCallback::assembleMatrix(const ID & matrix_id) { + if (matrix_id == "K") { + + } +} + +/* -------------------------------------------------------------------------- */ +void ContactSolverCallback::assembleLumpedMatrix(const ID & matrix_id) { + + +} + + +/* -------------------------------------------------------------------------- */ +void ContactSolverCallback::assembleResidual() { + /* ------------------------------------------------------------------------ */ + // computes the internal forces + solid.assembleInternalForces(); + solid.getDOFManager().assembleToResidual("displacement", + solid.getExternalForce(), 1); + solid.getDOFManager().assembleToResidual("displacement", + solid.getInternalForce(), 1); + +} + + +/* -------------------------------------------------------------------------- */ +MatrixType ContactSolverCallback::getMatrixType(const ID & matrix_id) { + + return _symmetric; +} + + +/* -------------------------------------------------------------------------- */ + int main(int argc, char *argv[]) { initialize("material_implicit.dat", argc, argv); const UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); mesh.read("hertz_implicit_2d.msh"); SolidMechanicsModel solid(mesh); solid.initFull(_analysis_method = _static); solid.setBaseName("static"); solid.addDumpFieldVector("displacement"); solid.addDumpField("blocked_dofs"); solid.addDumpField("external_force"); solid.addDumpField("internal_force"); solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top_body"); solid.applyBC(BC::Dirichlet::IncrementValue(-0.001, _y), "top_body"); solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bottom"); solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "bottom"); auto & solver = solid.getNonLinearSolver(); solver.set("max_iterations", 1000); solver.set("threshold", 1e-8); solver.set("convergence_type", _scc_residual); - auto & solid_stiffness = const_cast(solid.getDOFManager().getNewMatrix("K", _symmetric)); - solid.assembleInternalForces(); + //solid.assembleInternalForces(); solid.assembleStiffnessMatrix(); - //solid.solveStep(); - //solid.dump(); - + auto current_position = solid.getCurrentPosition(); ContactMechanicsModel contact(mesh, current_position); contact.initFull(_analysis_method = _implicit_contact); contact.setBaseNameToDumper("contact_mechanics", "contact"); contact.addDumpFieldVectorToDumper("contact_mechanics", "contact_force"); contact.addDumpFieldVectorToDumper("contact_mechanics", "external_force"); contact.addDumpFieldToDumper("contact_mechanics", "gaps"); contact.addDumpFieldToDumper("contact_mechanics", "areas"); contact.search(); contact.assembleInternalForces(); contact.assembleStiffnessMatrix(); contact.dump("paraview_all"); contact.dump("contact_mechanics"); Array & contact_force = contact.getInternalForce(); Array & external_force = solid.getExternalForce(); Array & blocked_dofs = solid.getBlockedDOFs(); for (auto && values: zip(make_view(external_force), make_view(contact_force), make_view(blocked_dofs)) ) { auto & ext_f = std::get<0>(values); auto & cont_f = std::get<1>(values); auto & blocked = std::get<2>(values); if (!blocked) { ext_f = cont_f; } - } + ContactSolverCallback callback(solid, contact); + auto & contact_stiffness = const_cast(contact.getDOFManager().getMatrix("K")); //auto & solid_stiffness = // const_cast(solid.getDOFManager().getNewMatrix("K", _symmetric)); solid_stiffness.add(contact_stiffness); - solid.solveStep(); + solid.solveStep(callback); contact.dump("paraview_all"); finalize(); return EXIT_SUCCESS; } diff --git a/src/common/aka_grid_dynamic.hh b/src/common/aka_grid_dynamic.hh index 9d7550fc4..1245cf799 100644 --- a/src/common/aka_grid_dynamic.hh +++ b/src/common/aka_grid_dynamic.hh @@ -1,513 +1,513 @@ /** * @file aka_grid_dynamic.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Wed Nov 08 2017 * * @brief Grid that is auto balanced * * @section LICENSE * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_common.hh" #include "aka_types.hh" #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_GRID_DYNAMIC_HH__ #define __AKANTU_AKA_GRID_DYNAMIC_HH__ namespace akantu { class Mesh; template class SpatialGrid { public: explicit SpatialGrid(UInt dimension) : dimension(dimension), spacing(dimension), center(dimension), lower(dimension), upper(dimension), empty_cell() {} SpatialGrid(UInt dimension, const Vector & spacing, const Vector & center) : dimension(dimension), spacing(spacing), center(center), lower(dimension), upper(dimension), empty_cell() { for (UInt i = 0; i < dimension; ++i) { lower(i) = std::numeric_limits::max(); upper(i) = -std::numeric_limits::max(); } } virtual ~SpatialGrid() = default; class neighbor_cells_iterator; class cells_iterator; class CellID { public: CellID() : ids() {} explicit CellID(UInt dimention) : ids(dimention) {} void setID(UInt dir, Int id) { ids(dir) = id; } Int getID(UInt dir) const { return ids(dir); } bool operator<(const CellID & id) const { return std::lexicographical_compare( ids.storage(), ids.storage() + ids.size(), id.ids.storage(), id.ids.storage() + id.ids.size()); } bool operator==(const CellID & id) const { return std::equal(ids.storage(), ids.storage() + ids.size(), id.ids.storage()); } bool operator!=(const CellID & id) const { return !(operator==(id)); } class neighbor_cells_iterator : private std::iterator { public: neighbor_cells_iterator(const CellID & cell_id, bool end) : cell_id(cell_id), position(cell_id.ids.size(), end ? 1 : -1) { this->updateIt(); if (end) this->it++; } neighbor_cells_iterator & operator++() { UInt i = 0; for (; i < position.size() && position(i) == 1; ++i) ; if (i == position.size()) { ++it; return *this; } for (UInt j = 0; j < i; ++j) position(j) = -1; position(i)++; updateIt(); return *this; } neighbor_cells_iterator operator++(int) { neighbor_cells_iterator tmp(*this); operator++(); return tmp; }; bool operator==(const neighbor_cells_iterator & rhs) const { return cell_id == rhs.cell_id && it == rhs.it; }; bool operator!=(const neighbor_cells_iterator & rhs) const { return !operator==(rhs); }; CellID operator*() const { CellID cur_cell_id(cell_id); cur_cell_id.ids += position; return cur_cell_id; }; private: void updateIt() { it = 0; for (UInt i = 0; i < position.size(); ++i) it = it * 3 + (position(i) + 1); } private: /// central cell id const CellID & cell_id; // number representing the current neighbor in base 3; UInt it; // current cell shift Vector position; }; class Neighbors { public: explicit Neighbors(const CellID & cell_id) : cell_id(cell_id) {} decltype(auto) begin() { return neighbor_cells_iterator(cell_id, false); } decltype(auto) end() { return neighbor_cells_iterator(cell_id, true); } private: const CellID & cell_id; }; decltype(auto) neighbors() { return Neighbors(*this); } private: friend class cells_iterator; Vector ids; }; /* ------------------------------------------------------------------------ */ class Cell { public: using iterator = typename std::vector::iterator; using const_iterator = typename std::vector::const_iterator; Cell() : id(), data() {} explicit Cell(const CellID & cell_id) : id(cell_id), data() {} bool operator==(const Cell & cell) const { return id == cell.id; } bool operator!=(const Cell & cell) const { return id != cell.id; } Cell & add(const T & d) { data.push_back(d); return *this; } iterator begin() { return data.begin(); } const_iterator begin() const { return data.begin(); } iterator end() { return data.end(); } const_iterator end() const { return data.end(); } private: CellID id; std::vector data; }; private: using cells_container = std::map; public: const Cell & getCell(const CellID & cell_id) const { auto it = cells.find(cell_id); if (it != cells.end()) return it->second; return empty_cell; } decltype(auto) beginCell(const CellID & cell_id) { auto it = cells.find(cell_id); if (it != cells.end()) return it->second.begin(); return empty_cell.begin(); } decltype(auto) endCell(const CellID & cell_id) { auto it = cells.find(cell_id); if (it != cells.end()) return it->second.end(); return empty_cell.end(); } decltype(auto) beginCell(const CellID & cell_id) const { auto it = cells.find(cell_id); if (it != cells.end()) return it->second.begin(); return empty_cell.begin(); } decltype(auto) endCell(const CellID & cell_id) const { auto it = cells.find(cell_id); if (it != cells.end()) return it->second.end(); return empty_cell.end(); } /* ------------------------------------------------------------------------ */ class cells_iterator : private std::iterator { public: explicit cells_iterator(typename std::map::const_iterator it) : it(it) {} cells_iterator & operator++() { this->it++; return *this; } cells_iterator operator++(int) { cells_iterator tmp(*this); operator++(); return tmp; }; bool operator==(const cells_iterator & rhs) const { return it == rhs.it; }; bool operator!=(const cells_iterator & rhs) const { return !operator==(rhs); }; CellID operator*() const { CellID cur_cell_id(this->it->first); return cur_cell_id; }; private: /// map iterator typename std::map::const_iterator it; }; public: template Cell & insert(const T & d, const vector_type & position) { auto && cell_id = getCellID(position); auto && it = cells.find(cell_id); if (it == cells.end()) { Cell cell(cell_id); auto & tmp = (cells[cell_id] = cell).add(d); for (UInt i = 0; i < dimension; ++i) { Real posl = center(i) + cell_id.getID(i) * spacing(i); Real posu = posl + spacing(i); - if (posl < lower(i)) + if (posl <= lower(i)) lower(i) = posl; if (posu > upper(i)) upper(i) = posu; } return tmp; } else { return it->second.add(d); } } /* ------------------------------------------------------------------------ */ inline decltype(auto) begin() const { auto begin = this->cells.begin(); return cells_iterator(begin); } inline decltype(auto) end() const { auto end = this->cells.end(); return cells_iterator(end); } template CellID getCellID(const vector_type & position) const { CellID cell_id(dimension); for (UInt i = 0; i < dimension; ++i) { cell_id.setID(i, getCellID(position(i), i)); } return cell_id; } void printself(std::ostream & stream, int indent = 0) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; std::streamsize prec = stream.precision(); std::ios_base::fmtflags ff = stream.flags(); stream.setf(std::ios_base::showbase); stream.precision(5); stream << space << "SpatialGrid<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; stream << space << " + dimension : " << this->dimension << std::endl; stream << space << " + lower bounds : {"; for (UInt i = 0; i < lower.size(); ++i) { if (i != 0) stream << ", "; stream << lower(i); }; stream << "}" << std::endl; stream << space << " + upper bounds : {"; for (UInt i = 0; i < upper.size(); ++i) { if (i != 0) stream << ", "; stream << upper(i); }; stream << "}" << std::endl; stream << space << " + spacing : {"; for (UInt i = 0; i < spacing.size(); ++i) { if (i != 0) stream << ", "; stream << spacing(i); }; stream << "}" << std::endl; stream << space << " + center : {"; for (UInt i = 0; i < center.size(); ++i) { if (i != 0) stream << ", "; stream << center(i); }; stream << "}" << std::endl; stream << space << " + nb_cells : " << this->cells.size() << "/"; Vector dist(this->dimension); dist = upper; dist -= lower; for (UInt i = 0; i < this->dimension; ++i) { dist(i) /= spacing(i); } UInt nb_cells = std::ceil(dist(0)); for (UInt i = 1; i < this->dimension; ++i) { nb_cells *= std::ceil(dist(i)); } stream << nb_cells << std::endl; stream << space << "]" << std::endl; stream.precision(prec); stream.flags(ff); } void saveAsMesh(Mesh & mesh) const; private: /* -------------------------------------------------------------------------- */ inline UInt getCellID(Real position, UInt direction) const { AKANTU_DEBUG_ASSERT(direction < center.size(), "The direction asked (" << direction << ") is out of range " << center.size()); Real dist_center = position - center(direction); Int id = std::floor(dist_center / spacing(direction)); // if(dist_center < 0) id--; return id; } friend class GridSynchronizer; public: AKANTU_GET_MACRO(LowerBounds, lower, const Vector &); AKANTU_GET_MACRO(UpperBounds, upper, const Vector &); AKANTU_GET_MACRO(Spacing, spacing, const Vector &); protected: UInt dimension; cells_container cells; Vector spacing; Vector center; Vector lower; Vector upper; Cell empty_cell; }; /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const SpatialGrid & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "mesh.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template void SpatialGrid::saveAsMesh(Mesh & mesh) const { auto & nodes = const_cast &>(mesh.getNodes()); ElementType type = _not_defined; switch (dimension) { case 1: type = _segment_2; break; case 2: type = _quadrangle_4; break; case 3: type = _hexahedron_8; break; } mesh.addConnectivityType(type); auto & connectivity = const_cast &>(mesh.getConnectivity(type)); auto & uint_data = mesh.getDataPointer("tag_1", type); Vector pos(dimension); UInt global_id = 0; for (auto & cell_pair : cells) { UInt cur_node = nodes.size(); UInt cur_elem = connectivity.size(); const CellID & cell_id = cell_pair.first; for (UInt i = 0; i < dimension; ++i) pos(i) = center(i) + cell_id.getID(i) * spacing(i); nodes.push_back(pos); for (UInt i = 0; i < dimension; ++i) pos(i) += spacing(i); nodes.push_back(pos); connectivity.push_back(cur_node); switch (dimension) { case 1: connectivity(cur_elem, 1) = cur_node + 1; break; case 2: pos(0) -= spacing(0); nodes.push_back(pos); pos(0) += spacing(0); pos(1) -= spacing(1); nodes.push_back(pos); connectivity(cur_elem, 1) = cur_node + 3; connectivity(cur_elem, 2) = cur_node + 1; connectivity(cur_elem, 3) = cur_node + 2; break; case 3: pos(1) -= spacing(1); pos(2) -= spacing(2); nodes.push_back(pos); pos(1) += spacing(1); nodes.push_back(pos); pos(0) -= spacing(0); nodes.push_back(pos); pos(1) -= spacing(1); pos(2) += spacing(2); nodes.push_back(pos); pos(0) += spacing(0); nodes.push_back(pos); pos(0) -= spacing(0); pos(1) += spacing(1); nodes.push_back(pos); connectivity(cur_elem, 1) = cur_node + 2; connectivity(cur_elem, 2) = cur_node + 3; connectivity(cur_elem, 3) = cur_node + 4; connectivity(cur_elem, 4) = cur_node + 5; connectivity(cur_elem, 5) = cur_node + 6; connectivity(cur_elem, 6) = cur_node + 1; connectivity(cur_elem, 7) = cur_node + 7; break; } uint_data.push_back(global_id); ++global_id; } } } // namespace akantu #endif /* __AKANTU_AKA_GRID_DYNAMIC_HH__ */ diff --git a/src/common/aka_math_tmpl.hh b/src/common/aka_math_tmpl.hh index 7af0c4683..2a825a001 100644 --- a/src/common/aka_math_tmpl.hh +++ b/src/common/aka_math_tmpl.hh @@ -1,783 +1,784 @@ /** * @file aka_math_tmpl.hh * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Alejandro M. Aragón * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Mathilde Radiguet * @author Nicolas Richart * @author Leonardo Snozzi * @author Peter Spijker * @author Marco Vocialta * * @date creation: Wed Aug 04 2010 * @date last modification: Tue Feb 20 2018 * * @brief Implementation of the inline functions of the math toolkit * * @section LICENSE * * Copyright (©) 2010-2018 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 . * */ } // akantu #include #include #include #include "aka_blas_lapack.hh" namespace akantu { /* -------------------------------------------------------------------------- */ inline void Math::matrix_vector(UInt im, UInt in, Real * A, Real * x, Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y char tran_A = 'N'; int incx = 1; int incy = 1; double beta = 0.; int m = im; int n = in; aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy); #else memset(y, 0, im * sizeof(Real)); for (UInt i = 0; i < im; ++i) { for (UInt j = 0; j < in; ++j) { y[i] += A[i + j * im] * x[j]; } y[i] *= alpha; } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_vector(UInt im, UInt in, Real * A, Real * x, Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y char tran_A = 'T'; int incx = 1; int incy = 1; double beta = 0.; int m = im; int n = in; aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy); #else memset(y, 0, in * sizeof(Real)); for (UInt i = 0; i < im; ++i) { for (UInt j = 0; j < in; ++j) { y[j] += A[j * im + i] * x[i]; } y[i] *= alpha; } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrix(UInt im, UInt in, UInt ik, Real * A, Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'N'; char trans_b = 'N'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &k, &beta, C, &m); #else memset(C, 0, im * in * sizeof(Real)); for (UInt j = 0; j < in; ++j) { UInt _jb = j * ik; UInt _jc = j * im; for (UInt i = 0; i < im; ++i) { for (UInt l = 0; l < ik; ++l) { UInt _la = l * im; C[i + _jc] += A[i + _la] * B[l + _jb]; } C[i + _jc] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrix(UInt im, UInt in, UInt ik, Real * A, Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'T'; char trans_b = 'N'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &k, &beta, C, &m); #else memset(C, 0, im * in * sizeof(Real)); for (UInt j = 0; j < in; ++j) { UInt _jc = j * im; UInt _jb = j * ik; for (UInt i = 0; i < im; ++i) { UInt _ia = i * ik; for (UInt l = 0; l < ik; ++l) { C[i + _jc] += A[l + _ia] * B[l + _jb]; } C[i + _jc] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrixt(UInt im, UInt in, UInt ik, Real * A, Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'N'; char trans_b = 'T'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &n, &beta, C, &m); #else memset(C, 0, im * in * sizeof(Real)); for (UInt j = 0; j < in; ++j) { UInt _jc = j * im; for (UInt i = 0; i < im; ++i) { for (UInt l = 0; l < ik; ++l) { UInt _la = l * im; UInt _lb = l * in; C[i + _jc] += A[i + _la] * B[j + _lb]; } C[i + _jc] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrixt(UInt im, UInt in, UInt ik, Real * A, Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'T'; char trans_b = 'T'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &n, &beta, C, &m); #else memset(C, 0, im * in * sizeof(Real)); for (UInt j = 0; j < in; ++j) { UInt _jc = j * im; for (UInt i = 0; i < im; ++i) { UInt _ia = i * ik; for (UInt l = 0; l < ik; ++l) { UInt _lb = l * in; C[i + _jc] += A[l + _ia] * B[j + _lb]; } C[i + _jc] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::aXplusY(UInt n, Real alpha, Real * x, Real * y) { #ifdef AKANTU_USE_BLAS /// y := alpha x + y int incx = 1, incy = 1; aka_axpy(&n, &alpha, x, &incx, y, &incy); #else for (UInt i = 0; i < n; ++i) *(y++) += alpha * *(x++); #endif } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot(Real * v1, Real * v2, UInt in) { #ifdef AKANTU_USE_BLAS /// d := v1 . v2 int incx = 1, incy = 1, n = in; Real d = aka_dot(&n, v1, &incx, v2, &incy); #else Real d = 0; for (UInt i = 0; i < in; ++i) { d += v1[i] * v2[i]; } #endif return d; } /* -------------------------------------------------------------------------- */ template inline void Math::matMul(UInt m, UInt n, UInt k, Real alpha, Real * A, Real * B, __attribute__((unused)) Real beta, Real * C) { if (tr_A) { if (tr_B) matrixt_matrixt(m, n, k, A, B, C, alpha); else matrixt_matrix(m, n, k, A, B, C, alpha); } else { if (tr_B) matrix_matrixt(m, n, k, A, B, C, alpha); else matrix_matrix(m, n, k, A, B, C, alpha); } } /* -------------------------------------------------------------------------- */ template inline void Math::matVectMul(UInt m, UInt n, Real alpha, Real * A, Real * x, __attribute__((unused)) Real beta, Real * y) { if (tr_A) { matrixt_vector(m, n, A, x, y, alpha); } else { matrix_vector(m, n, A, x, y, alpha); } } /* -------------------------------------------------------------------------- */ template inline void Math::matrixEig(UInt n, T * A, T * d, T * V) { // Matrix A is row major, so the lapack function in fortran will process // A^t. Asking for the left eigenvectors of A^t will give the transposed right // eigenvectors of A so in the C++ code the right eigenvectors. char jobvr, jobvl; if (V != nullptr) jobvr = 'V'; // compute left eigenvectors else jobvr = 'N'; // compute left eigenvectors jobvl = 'N'; // compute right eigenvectors auto * di = new T[n]; // imaginary part of the eigenvalues int info; int N = n; T wkopt; int lwork = -1; // query and allocate the optimal workspace aka_geev(&jobvl, &jobvr, &N, A, &N, d, di, nullptr, &N, V, &N, &wkopt, &lwork, &info); lwork = int(wkopt); auto * work = new T[lwork]; // solve the eigenproblem aka_geev(&jobvl, &jobvr, &N, A, &N, d, di, nullptr, &N, V, &N, work, &lwork, &info); AKANTU_DEBUG_ASSERT( info == 0, "Problem computing eigenvalues/vectors. DGEEV exited with the value " << info); delete[] work; delete[] di; // I hope for you that there was no complex eigenvalues !!! } /* -------------------------------------------------------------------------- */ inline void Math::matrix22_eigenvalues(Real * A, Real * Adiag) { /// d = determinant of Matrix A Real d = det2(A); /// b = trace of Matrix A Real b = A[0] + A[3]; Real c = sqrt(b * b - 4 * d); Adiag[0] = .5 * (b + c); Adiag[1] = .5 * (b - c); } /* -------------------------------------------------------------------------- */ inline void Math::matrix33_eigenvalues(Real * A, Real * Adiag) { matrixEig(3, A, Adiag); } /* -------------------------------------------------------------------------- */ template inline void Math::eigenvalues(Real * A, Real * d) { if (dim == 1) { d[0] = A[0]; } else if (dim == 2) { matrix22_eigenvalues(A, d); } // else if(dim == 3) { matrix33_eigenvalues(A, d); } else matrixEig(dim, A, d); } /* -------------------------------------------------------------------------- */ inline Real Math::det2(const Real * mat) { return mat[0] * mat[3] - mat[1] * mat[2]; } /* -------------------------------------------------------------------------- */ inline Real Math::det3(const Real * mat) { return mat[0] * (mat[4] * mat[8] - mat[7] * mat[5]) - mat[3] * (mat[1] * mat[8] - mat[7] * mat[2]) + mat[6] * (mat[1] * mat[5] - mat[4] * mat[2]); } /* -------------------------------------------------------------------------- */ template inline Real Math::det(const Real * mat) { if (n == 1) return *mat; else if (n == 2) return det2(mat); else if (n == 3) return det3(mat); else return det(n, mat); } /* -------------------------------------------------------------------------- */ template inline T Math::det(UInt n, const T * A) { int N = n; int info; auto * ipiv = new int[N + 1]; auto * LU = new T[N * N]; std::copy(A, A + N * N, LU); // LU factorization of A aka_getrf(&N, &N, LU, &N, ipiv, &info); if (info > 0) { AKANTU_ERROR("Singular matrix - cannot factorize it (info: " << info << " )"); } // det(A) = det(L) * det(U) = 1 * det(U) = product_i U_{ii} T det = 1.; for (int i = 0; i < N; ++i) det *= (2 * (ipiv[i] == i) - 1) * LU[i * n + i]; delete[] ipiv; delete[] LU; return det; } /* -------------------------------------------------------------------------- */ inline void Math::normal2(const Real * vec, Real * normal) { normal[0] = vec[1]; normal[1] = -vec[0]; Math::normalize2(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normal3(const Real * vec1, const Real * vec2, Real * normal) { Math::vectorProduct3(vec1, vec2, normal); Math::normalize3(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normalize2(Real * vec) { Real norm = Math::norm2(vec); vec[0] /= norm; vec[1] /= norm; } /* -------------------------------------------------------------------------- */ inline void Math::normalize3(Real * vec) { Real norm = Math::norm3(vec); vec[0] /= norm; vec[1] /= norm; vec[2] /= norm; } /* -------------------------------------------------------------------------- */ inline Real Math::norm2(const Real * vec) { return sqrt(vec[0] * vec[0] + vec[1] * vec[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::norm3(const Real * vec) { return sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); } /* -------------------------------------------------------------------------- */ inline Real Math::norm(UInt n, const Real * vec) { Real norm = 0.; for (UInt i = 0; i < n; ++i) { norm += vec[i] * vec[i]; } return sqrt(norm); } /* -------------------------------------------------------------------------- */ inline void Math::inv2(const Real * mat, Real * inv) { Real det_mat = det2(mat); inv[0] = mat[3] / det_mat; inv[1] = -mat[1] / det_mat; inv[2] = -mat[2] / det_mat; inv[3] = mat[0] / det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::inv3(const Real * mat, Real * inv) { Real det_mat = det3(mat); inv[0] = (mat[4] * mat[8] - mat[7] * mat[5]) / det_mat; inv[1] = (mat[2] * mat[7] - mat[8] * mat[1]) / det_mat; inv[2] = (mat[1] * mat[5] - mat[4] * mat[2]) / det_mat; inv[3] = (mat[5] * mat[6] - mat[8] * mat[3]) / det_mat; inv[4] = (mat[0] * mat[8] - mat[6] * mat[2]) / det_mat; inv[5] = (mat[2] * mat[3] - mat[5] * mat[0]) / det_mat; inv[6] = (mat[3] * mat[7] - mat[6] * mat[4]) / det_mat; inv[7] = (mat[1] * mat[6] - mat[7] * mat[0]) / det_mat; inv[8] = (mat[0] * mat[4] - mat[3] * mat[1]) / det_mat; } /* -------------------------------------------------------------------------- */ template inline void Math::inv(const Real * A, Real * Ainv) { if (n == 1) *Ainv = 1. / *A; else if (n == 2) inv2(A, Ainv); else if (n == 3) inv3(A, Ainv); else inv(n, A, Ainv); } /* -------------------------------------------------------------------------- */ template inline void Math::inv(UInt n, const T * A, T * invA) { int N = n; int info; auto * ipiv = new int[N + 1]; int lwork = N * N; auto * work = new T[lwork]; std::copy(A, A + n * n, invA); aka_getrf(&N, &N, invA, &N, ipiv, &info); if (info > 0) { AKANTU_ERROR("Singular matrix - cannot factorize it (info: " << info << " )"); } aka_getri(&N, invA, &N, ipiv, work, &lwork, &info); if (info != 0) { AKANTU_ERROR("Cannot invert the matrix (info: " << info << " )"); } delete[] ipiv; delete[] work; } /* -------------------------------------------------------------------------- */ template inline void Math::solve(UInt n, const T * A, T * x, const T * b) { int N = n; int info; auto * ipiv = new int[N]; auto * lu_A = new T[N * N]; std::copy(A, A + N * N, lu_A); aka_getrf(&N, &N, lu_A, &N, ipiv, &info); if (info > 0) { AKANTU_ERROR("Singular matrix - cannot factorize it (info: " << info << " )"); } char trans = 'N'; int nrhs = 1; std::copy(b, b + N, x); aka_getrs(&trans, &N, &nrhs, lu_A, &N, ipiv, x, &N, &info); if (info != 0) { AKANTU_ERROR("Cannot solve the system (info: " << info << " )"); } delete[] ipiv; delete[] lu_A; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ inline Real Math::matrixDoubleDot22(Real * A, Real * B) { Real d; d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3]; return d; } /* -------------------------------------------------------------------------- */ inline Real Math::matrixDoubleDot33(Real * A, Real * B) { Real d; d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3] + A[4] * B[4] + A[5] * B[5] + A[6] * B[6] + A[7] * B[7] + A[8] * B[8]; return d; } /* -------------------------------------------------------------------------- */ inline Real Math::matrixDoubleDot(UInt n, Real * A, Real * B) { Real d = 0.; for (UInt i = 0; i < n; ++i) { for (UInt j = 0; j < n; ++j) { d += A[i * n + j] * B[i * n + j]; } } return d; } /* -------------------------------------------------------------------------- */ inline void Math::vectorProduct3(const Real * v1, const Real * v2, Real * res) { res[0] = v1[1] * v2[2] - v1[2] * v2[1]; res[1] = v1[2] * v2[0] - v1[0] * v2[2]; res[2] = v1[0] * v2[1] - v1[1] * v2[0]; } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot2(const Real * v1, const Real * v2) { return (v1[0] * v2[0] + v1[1] * v2[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot3(const Real * v1, const Real * v2) { return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_2d(const Real * x, const Real * y) { return sqrt((y[0] - x[0]) * (y[0] - x[0]) + (y[1] - x[1]) * (y[1] - x[1])); } /* -------------------------------------------------------------------------- */ inline Real Math::triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3) { /** * @f{eqnarray*}{ * r &=& A / s \\ * A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)} \\ * s &=& \frac{a + b + c}{2} * @f} */ Real a, b, c; + a = distance_2d(coord1, coord2); b = distance_2d(coord2, coord3); c = distance_2d(coord1, coord3); Real s; s = (a + b + c) * 0.5; return sqrt((s - a) * (s - b) * (s - c) / s); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_3d(const Real * x, const Real * y) { return sqrt((y[0] - x[0]) * (y[0] - x[0]) + (y[1] - x[1]) * (y[1] - x[1]) + (y[2] - x[2]) * (y[2] - x[2])); } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real xx[9], vol; xx[0] = coord2[0]; xx[1] = coord2[1]; xx[2] = coord2[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol = det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol -= det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol += det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord3[0]; xx[7] = coord3[1]; xx[8] = coord3[2]; vol -= det3(xx); vol /= 6; return vol; } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real l12, l13, l14, l23, l24, l34; l12 = distance_3d(coord1, coord2); l13 = distance_3d(coord1, coord3); l14 = distance_3d(coord1, coord4); l23 = distance_3d(coord2, coord3); l24 = distance_3d(coord2, coord4); l34 = distance_3d(coord3, coord4); Real s1, s2, s3, s4; s1 = (l12 + l23 + l13) * 0.5; s1 = sqrt(s1 * (s1 - l12) * (s1 - l23) * (s1 - l13)); s2 = (l12 + l24 + l14) * 0.5; s2 = sqrt(s2 * (s2 - l12) * (s2 - l24) * (s2 - l14)); s3 = (l23 + l34 + l24) * 0.5; s3 = sqrt(s3 * (s3 - l23) * (s3 - l34) * (s3 - l24)); s4 = (l13 + l34 + l14) * 0.5; s4 = sqrt(s4 * (s4 - l13) * (s4 - l34) * (s4 - l14)); Real volume = Math::tetrahedron_volume(coord1, coord2, coord3, coord4); return 3 * volume / (s1 + s2 + s3 + s4); } /* -------------------------------------------------------------------------- */ inline void Math::barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter) { memset(barycenter, 0, spatial_dimension * sizeof(Real)); for (UInt n = 0; n < nb_points; ++n) { UInt offset = n * spatial_dimension; for (UInt i = 0; i < spatial_dimension; ++i) { barycenter[i] += coord[offset + i] / (Real)nb_points; } } } /* -------------------------------------------------------------------------- */ inline void Math::vector_2d(const Real * x, const Real * y, Real * res) { res[0] = y[0] - x[0]; res[1] = y[1] - x[1]; } /* -------------------------------------------------------------------------- */ inline void Math::vector_3d(const Real * x, const Real * y, Real * res) { res[0] = y[0] - x[0]; res[1] = y[1] - x[1]; res[2] = y[2] - x[2]; } /* -------------------------------------------------------------------------- */ /// Combined absolute and relative tolerance test proposed in /// Real-time collision detection by C. Ericson (2004) inline bool Math::are_float_equal(const Real x, const Real y) { Real abs_max = std::max(std::abs(x), std::abs(y)); abs_max = std::max(abs_max, Real(1.)); return std::abs(x - y) <= (tolerance * abs_max); } /* -------------------------------------------------------------------------- */ inline bool Math::isnan(Real x) { #if defined(__INTEL_COMPILER) #pragma warning(push) #pragma warning(disable : 1572) #endif // defined(__INTEL_COMPILER) // x = x return false means x = quiet_NaN return !(x == x); #if defined(__INTEL_COMPILER) #pragma warning(pop) #endif // defined(__INTEL_COMPILER) } /* -------------------------------------------------------------------------- */ inline bool Math::are_vector_equal(UInt n, Real * x, Real * y) { bool test = true; for (UInt i = 0; i < n; ++i) { test &= are_float_equal(x[i], y[i]); } return test; } /* -------------------------------------------------------------------------- */ inline bool Math::intersects(Real x_min, Real x_max, Real y_min, Real y_max) { return not ((x_max < y_min) or (x_min > y_max)); } /* -------------------------------------------------------------------------- */ inline bool Math::is_in_range(Real a, Real x_min, Real x_max) { return ((a >= x_min) and (a <= x_max)); } /* -------------------------------------------------------------------------- */ template inline T Math::pow(T x) { return (pow

(x) * x); } template <> inline UInt Math::pow<0, UInt>(__attribute__((unused)) UInt x) { return (1); } template <> inline Real Math::pow<0, Real>(__attribute__((unused)) Real x) { return (1.); } /* -------------------------------------------------------------------------- */ template Real Math::NewtonRaphson::solve(const Functor & funct, Real x_0) { Real x = x_0; Real f_x = funct.f(x); UInt iter = 0; while (std::abs(f_x) > this->tolerance && iter < this->max_iteration) { x -= f_x / funct.f_prime(x); f_x = funct.f(x); iter++; } AKANTU_DEBUG_ASSERT(iter < this->max_iteration, "Newton Raphson (" << funct.name << ") solve did not converge in " << this->max_iteration << " iterations (tolerance: " << this->tolerance << ")"); return x; } diff --git a/src/model/model_solver.cc b/src/model/#model_solver.cc# similarity index 98% copy from src/model/model_solver.cc copy to src/model/#model_solver.cc# index 1891e7fd3..ce889b9f9 100644 --- a/src/model/model_solver.cc +++ b/src/model/#model_solver.cc# @@ -1,365 +1,371 @@ /** * @file model_solver.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of ModelSolver * * @section LICENSE * * Copyright (©) 2015-2018 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 "model_solver.hh" #include "dof_manager.hh" #include "dof_manager_default.hh" #include "mesh.hh" #include "non_linear_solver.hh" #include "time_step_solver.hh" #if defined(AKANTU_USE_PETSC) #include "dof_manager_petsc.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template static T getOptionToType(const std::string & opt_str) { std::stringstream sstr(opt_str); T opt; sstr >> opt; return opt; } /* -------------------------------------------------------------------------- */ ModelSolver::ModelSolver(Mesh & mesh, const ModelType & type, const ID & id, UInt memory_id) : Parsable(ParserType::_model, id), SolverCallback(), model_type(type), parent_id(id), parent_memory_id(memory_id), mesh(mesh), dof_manager(nullptr), default_solver_id("") {} /* -------------------------------------------------------------------------- */ ModelSolver::~ModelSolver() = default; /* -------------------------------------------------------------------------- */ std::tuple ModelSolver::getParserSection() { auto sub_sections = getStaticParser().getSubSections(ParserType::_model); auto it = std::find_if( sub_sections.begin(), sub_sections.end(), [&](auto && section) { ModelType type = getOptionToType(section.getName()); // default id should be the model type if not defined std::string name = section.getParameter("name", this->parent_id); return type == model_type and name == this->parent_id; }); if (it == sub_sections.end()) return std::make_tuple(ParserSection(), true); return std::make_tuple(*it, false); } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager() { // default without external solver activated at compilation same as mumps that // is the historical solver but with only the lumped solver ID solver_type = "default"; #if defined(AKANTU_USE_MUMPS) solver_type = "default"; #elif defined(AKANTU_USE_PETSC) solver_type = "petsc"; #endif ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { solver_type = section.getOption(solver_type); this->initDOFManager(section, solver_type); } else { this->initDOFManager(solver_type); } } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager(const ID & solver_type) { try { this->dof_manager = DOFManagerFactory::getInstance().allocate( solver_type, mesh, this->parent_id + ":dof_manager" + solver_type, this->parent_memory_id); } catch (...) { AKANTU_EXCEPTION( "To use the solver " << solver_type << " you will have to code it. This is an unknown solver type."); } this->setDOFManager(*this->dof_manager); } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager(const ParserSection & section, const ID & solver_type) { this->initDOFManager(solver_type); auto sub_sections = section.getSubSections(ParserType::_time_step_solver); // parsing the time step solvers for (auto && section : sub_sections) { ID type = section.getName(); ID solver_id = section.getParameter("name", type); auto tss_type = getOptionToType(type); auto tss_options = this->getDefaultSolverOptions(tss_type); auto sub_solvers_sect = section.getSubSections(ParserType::_non_linear_solver); auto nb_non_linear_solver_section = section.getNbSubSections(ParserType::_non_linear_solver); auto nls_type = tss_options.non_linear_solver_type; if (nb_non_linear_solver_section == 1) { auto && nls_section = *(sub_solvers_sect.first); nls_type = getOptionToType(nls_section.getName()); } else if (nb_non_linear_solver_section > 0) { AKANTU_EXCEPTION("More than one non linear solver are provided for the " "time step solver " << solver_id); } this->getNewSolver(solver_id, tss_type, nls_type); if (nb_non_linear_solver_section == 1) { const auto & nls_section = *(sub_solvers_sect.first); this->dof_manager->getNonLinearSolver(solver_id).parseSection( nls_section); } auto sub_integrator_sections = section.getSubSections(ParserType::_integration_scheme); for (auto && is_section : sub_integrator_sections) { const auto & dof_type_str = is_section.getName(); ID dof_id; try { ID tmp = is_section.getParameter("name"); dof_id = tmp; } catch (...) { AKANTU_EXCEPTION("No degree of freedom name specified for the " "integration scheme of type " << dof_type_str); } auto it_type = getOptionToType(dof_type_str); IntegrationScheme::SolutionType s_type = is_section.getParameter( "solution_type", tss_options.solution_type[dof_id]); this->setIntegrationScheme(solver_id, dof_id, it_type, s_type); } for (auto & is_type : tss_options.integration_scheme_type) { if (!this->hasIntegrationScheme(solver_id, is_type.first)) { this->setIntegrationScheme(solver_id, is_type.first, is_type.second, tss_options.solution_type[is_type.first]); } } } if (section.hasParameter("default_solver")) { ID default_solver = section.getParameter("default_solver"); if (this->hasSolver(default_solver)) { this->setDefaultSolver(default_solver); } else AKANTU_EXCEPTION( "The solver \"" << default_solver << "\" was not created, it cannot be set as default solver"); } } /* -------------------------------------------------------------------------- */ TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) { ID tmp_solver_id = solver_id; if (tmp_solver_id == "") tmp_solver_id = this->default_solver_id; TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id); return tss; } /* -------------------------------------------------------------------------- */ const TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) const { ID tmp_solver_id = solver_id; if (solver_id == "") tmp_solver_id = this->default_solver_id; const TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id); return tss; } /* -------------------------------------------------------------------------- */ TimeStepSolver & ModelSolver::getTimeStepSolver(const ID & solver_id) { return this->getSolver(solver_id); } /* -------------------------------------------------------------------------- */ const TimeStepSolver & ModelSolver::getTimeStepSolver(const ID & solver_id) const { return this->getSolver(solver_id); } /* -------------------------------------------------------------------------- */ NonLinearSolver & ModelSolver::getNonLinearSolver(const ID & solver_id) { return this->getSolver(solver_id).getNonLinearSolver(); } /* -------------------------------------------------------------------------- */ const NonLinearSolver & ModelSolver::getNonLinearSolver(const ID & solver_id) const { return this->getSolver(solver_id).getNonLinearSolver(); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasSolver(const ID & solver_id) const { ID tmp_solver_id = solver_id; if (solver_id == "") tmp_solver_id = this->default_solver_id; if (not this->dof_manager) { AKANTU_EXCEPTION("No DOF manager was initialized"); } return this->dof_manager->hasTimeStepSolver(tmp_solver_id); } /* -------------------------------------------------------------------------- */ void ModelSolver::setDefaultSolver(const ID & solver_id) { AKANTU_DEBUG_ASSERT( this->hasSolver(solver_id), "Cannot set the default solver to a solver that does not exists"); this->default_solver_id = solver_id; } /* -------------------------------------------------------------------------- */ void ModelSolver::solveStep(const ID & solver_id) { + solveStep(*this, solver_id); +} + +/* -------------------------------------------------------------------------- */ +void ModelSolver::solveStep(SolverCallback & callback, + const ID & solver_id) { AKANTU_DEBUG_IN(); TimeStepSolver & tss = this->getSolver(solver_id); // make one non linear solve - tss.solveStep(*this); + tss.solveStep(callback); AKANTU_DEBUG_OUT(); } - + /* -------------------------------------------------------------------------- */ void ModelSolver::getNewSolver(const ID & solver_id, TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { if (this->default_solver_id == "") { this->default_solver_id = solver_id; } if (non_linear_solver_type == _nls_auto) { switch (time_step_solver_type) { case _tsst_dynamic: case _tsst_static: non_linear_solver_type = _nls_newton_raphson; break; case _tsst_dynamic_lumped: non_linear_solver_type = _nls_lumped; break; case _tsst_not_defined: AKANTU_EXCEPTION(time_step_solver_type << " is not a valid time step solver type"); break; } } this->initSolver(time_step_solver_type, non_linear_solver_type); NonLinearSolver & nls = this->dof_manager->getNewNonLinearSolver( solver_id, non_linear_solver_type); this->dof_manager->getNewTimeStepSolver(solver_id, time_step_solver_type, nls); } /* -------------------------------------------------------------------------- */ Real ModelSolver::getTimeStep(const ID & solver_id) const { const TimeStepSolver & tss = this->getSolver(solver_id); return tss.getTimeStep(); } /* -------------------------------------------------------------------------- */ void ModelSolver::setTimeStep(Real time_step, const ID & solver_id) { TimeStepSolver & tss = this->getSolver(solver_id); return tss.setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void ModelSolver::setIntegrationScheme( const ID & solver_id, const ID & dof_id, const IntegrationSchemeType & integration_scheme_type, IntegrationScheme::SolutionType solution_type) { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); tss.setIntegrationScheme(dof_id, integration_scheme_type, solution_type); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasDefaultSolver() const { return (this->default_solver_id != ""); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasIntegrationScheme(const ID & solver_id, const ID & dof_id) const { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); return tss.hasIntegrationScheme(dof_id); } /* -------------------------------------------------------------------------- */ void ModelSolver::predictor() {} /* -------------------------------------------------------------------------- */ void ModelSolver::corrector() {} /* -------------------------------------------------------------------------- */ TimeStepSolverType ModelSolver::getDefaultSolverType() const { return _tsst_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions ModelSolver::getDefaultSolverOptions(__attribute__((unused)) const TimeStepSolverType & type) const { ModelSolverOptions options; options.non_linear_solver_type = _nls_auto; return options; } } // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index 6f5ada151..ae96e37c3 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,564 +1,615 @@ /** * @file contact_detector.cc * * @author Mohit Pundir * * @date creation: Wed Sep 12 2018 * @date last modification: Fri Sep 21 2018 * * @brief Mother class for all detection algorithms * * @section LICENSE * * Copyright (©) 2010-2018 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 "contact_detector.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactDetector::ContactDetector(Mesh & mesh, const ID & id, UInt memory_id) : ContactDetector(mesh, mesh.getNodes(), id, memory_id) { } /* -------------------------------------------------------------------------- */ - ContactDetector::ContactDetector(Mesh & mesh, Array & positions, const ID & id, UInt memory_id) +ContactDetector::ContactDetector(Mesh & mesh, Array & positions, const ID & id, UInt memory_id) : Memory(id, memory_id), Parsable(ParserType::_contact_detector, id), mesh(mesh), positions(positions) { AKANTU_DEBUG_IN(); this->spatial_dimension = mesh.getSpatialDimension(); this->mesh.fillNodesToElements(this->spatial_dimension - 1); this->parseSection(); this->getMaximalDetectionDistance(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactDetector::parseSection() { const Parser & parser = getStaticParser(); const ParserSection & section = *(parser.getSubSections(ParserType::_contact_detector).first); auto type = section.getParameterValue("type"); if (type == "implicit") { this->detection_type = _implicit; } else if (type == "explicit"){ this->detection_type = _explicit; } else { AKANTU_ERROR("Unknown detection type : " << type); } surfaces[Surface::master] = section.getParameterValue("master_surface"); surfaces[Surface::slave ] = section.getParameterValue("slave_surface"); } /* -------------------------------------------------------------------------- */ void ContactDetector::getMaximalDetectionDistance() { AKANTU_DEBUG_IN(); - Real el_size; - Real max_el_size = std::numeric_limits::min(); + Real elem_size; + Real max_elem_size = std::numeric_limits::min(); + std::cerr << max_elem_size << std::endl; auto & master_group = mesh.getElementGroup(surfaces[Surface::master]); - for (auto & type: master_group.elementTypes(spatial_dimension - 1, _not_ghost)) { - + for (auto type: + master_group.elementTypes(spatial_dimension - 1, _not_ghost, _ek_regular)) { + + const auto & element_ids = master_group.getElements(type); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); + UInt nb_elements = element_ids.size(); + + Element elem; + elem.type = type; + for (auto el : element_ids) { + elem.element = el; + Matrix elem_coords(spatial_dimension, nb_nodes_per_element); + this->coordinatesOfElement(elem, elem_coords); - Array coord(0, nb_nodes_per_element * spatial_dimension); - FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, - _not_ghost); - auto el_coord = coord.begin(spatial_dimension, nb_nodes_per_element); - UInt nb_element = mesh.getNbElement(type); + elem_size = FEEngine::getElementInradius(elem_coords, type); - for (UInt el =0; el < nb_element; ++el, ++el_coord) { - el_size = FEEngine::getElementInradius(*el_coord, type); - max_el_size = std::max(max_el_size, el_size); + std::cerr << elem_coords << std::endl; + std::cerr << elem_size << std::endl; + + max_elem_size = std::max(max_elem_size, elem_size); } AKANTU_DEBUG_INFO("The maximum element size : " - << max_el_size ); + << max_elem_size ); + + std::cerr << max_elem_size << std::endl; } - this->max_dd = max_el_size; - this->max_bb = max_el_size; + this->max_dd = max_elem_size; + this->max_bb = max_elem_size; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactDetector::search(std::map & contact_map) { this->globalSearch(contact_map); } /* -------------------------------------------------------------------------- */ void ContactDetector::globalSearch(std::map & contact_map) { auto & master_list = mesh.getElementGroup(surfaces[Surface::master]).getNodeGroup().getNodes(); auto & slave_list = mesh.getElementGroup(surfaces[Surface::slave]).getNodeGroup().getNodes(); BBox bbox_master(spatial_dimension); this->constructBoundingBox(bbox_master, master_list); BBox bbox_slave(spatial_dimension); this->constructBoundingBox(bbox_slave, slave_list); auto && bbox_intersection = bbox_master.intersection(bbox_slave); + + std::cerr << bbox_master << std::endl; + std::cerr << bbox_slave << std::endl; + + std::cerr << bbox_intersection << std::endl; AKANTU_DEBUG_INFO( "Intersection BBox " << bbox_intersection ); Vector center(spatial_dimension); bbox_intersection.getCenter(center); Vector spacing(spatial_dimension); this->computeCellSpacing(spacing); auto & master_surface_list = mesh.getElementGroup(surfaces[Surface::master]).getNodeGroup().getNodes(); auto & slave_surface_list = mesh.getElementGroup(surfaces[Surface::slave]).getNodeGroup().getNodes(); SpatialGrid master_grid(spatial_dimension, spacing, center); this->constructGrid(master_grid, bbox_intersection, master_surface_list); SpatialGrid slave_grid(spatial_dimension, spacing, center); this->constructGrid(slave_grid, bbox_intersection, slave_surface_list); if (AKANTU_DEBUG_TEST(dblDump)) { Mesh mesh(spatial_dimension, "save"); master_grid.saveAsMesh(mesh); mesh.write("master_grid.msh"); } if (AKANTU_DEBUG_TEST(dblDump)) { Mesh mesh2(spatial_dimension, "save"); slave_grid.saveAsMesh(mesh2); mesh2.write("slave_grid.msh"); } AKANTU_DEBUG_INFO( "Grid Details " << master_grid ); // search slave grid nodes in contactelement array and if they exits // and still have orthogonal projection on its associated master // facetremove it from the spatial grid or do not consider it for // local search, maybe better option will be to have spatial grid of // type node info and one of the variable of node info should be // facet already exits // so contact eleemnts will be updated based on the above // consideration , this means only those contact elements will be // keep whose slave node is still in intersection bbox and still has // projection in its master facet // also if slave node is already exists in contact element and // orthogonal projection does not exits then search the associated // master facets with the current master facets within a given // radius , this is subjected to computational cost as searching // neighbbor cells can be more effective. this->localSearch(slave_grid, master_grid, contact_map); } /* -------------------------------------------------------------------------- */ void ContactDetector::localSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid, std::map & contact_map) { // local search // out of these array check each cell for closet node in that cell // and neighbouring cells find the actual orthogonally closet // check the projection of slave node on master facets connected to // the closet master node, if yes update the contact element with // slave node and master node and master surfaces connected to the // master node // these master surfaces will be needed later to update contact // elements Array slave_nodes; Array master_nodes; BBox bbox_master_grid(spatial_dimension); BBox bbox_slave_grid(spatial_dimension); - bbox_master_grid += master_grid.getUpperBounds(); - bbox_master_grid += master_grid.getLowerBounds(); + auto create_bbox = [&](auto & grid, auto & bbox) { + auto upper_bound = grid.getUpperBounds(); + auto lower_bound = grid.getLowerBounds(); + for (UInt s: arange(spatial_dimension)) { + lower_bound(s) -= this->max_bb; + upper_bound(s) += this->max_bb; + } - bbox_slave_grid += slave_grid.getUpperBounds(); - bbox_slave_grid += slave_grid.getLowerBounds(); + bbox += lower_bound; + bbox += upper_bound; + }; + create_bbox(master_grid, bbox_master_grid); + create_bbox(slave_grid, bbox_slave_grid); + auto && bbox_intersection = bbox_master_grid.intersection(bbox_slave_grid); + + std::cerr << bbox_master_grid << std::endl; + std::cerr << bbox_slave_grid << std::endl; + + std::cerr << bbox_intersection << std::endl; // find the closet master node for each slave node for (auto && cell_id : slave_grid) { AKANTU_DEBUG_INFO("Looping on next cell"); for (auto && q1: slave_grid.getCell(cell_id)) { - + + bool pair_exists = false; + Vector pos(spatial_dimension); for (UInt s: arange(spatial_dimension)) { pos(s) = this->positions(q1, s); } if (!bbox_intersection.contains(pos)) { continue; } - Real closet_distance = std::numeric_limits::max(); + Real closet_distance = this->max_dd * 0.5; UInt closet_master_node; // loop over all the neighboring cells of the current cells for (auto && neighbor_cell : cell_id.neighbors()) { // loop over the data of neighboring cells from master grid for (auto && q2 : master_grid.getCell(neighbor_cell)) { AKANTU_DEBUG_INFO("Looping on neighbor cell in master"); Vector pos2(spatial_dimension); for (UInt s: arange(spatial_dimension)) { pos2(s) = this->positions(q2, s); } Real distance = pos.distance(pos2); if (distance <= closet_distance) { closet_master_node = q2; closet_distance = distance; + pair_exists = true; } } } - - slave_nodes.push_back(q1); - master_nodes.push_back(closet_master_node); + + if (pair_exists) { + slave_nodes.push_back(q1); + master_nodes.push_back(closet_master_node); + std::cerr << q1 << " ---- " << closet_master_node << "---" << closet_distance << std::endl; + } } } + + std::cerr << "Number of Slave nodes = " << slave_nodes.size() << std::endl; for (auto && values : zip(slave_nodes, master_nodes)) { const auto & slave_node = std::get<0>(values); const auto & master_node = std::get<1>(values); Array elements; this->mesh.getAssociatedElements(master_node, elements); auto normals = std::make_unique>(elements.size(), spatial_dimension, "normals"); auto gaps = std::make_unique>(elements.size(), 1, "gaps"); auto natural_projections = std::make_unique>(elements.size(), spatial_dimension - 1, "projections"); - + auto status = std::make_unique>(elements.size(), + 1, "status"); this->computeOrthogonalProjection(slave_node, elements, - *normals, *gaps, *natural_projections); - - auto minimum = std::min_element( gaps->begin(), gaps->end()); - auto index = std::distance( gaps->begin(), minimum); - + *normals, *gaps, *natural_projections, *status); + + + UInt index; + Real minimum_gap = std::numeric_limits::max(); + bool to_consider = false; + for (UInt i : arange(gaps->size())) { + if (!(*status)[i]) + continue; + + if ((*gaps)[i] <= minimum_gap) { + minimum_gap = (*gaps)[i]; + index = i; + to_consider = true; + } + } + + if (!to_consider) { + continue; + } + Vector master_conn = this->mesh.getConnectivity(elements[index]); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave_node; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i-1]; } - + contact_map[slave_node] = ContactElement(elements[index]); contact_map[slave_node].gap = (*gaps)[index]; contact_map[slave_node].normal = Vector(normals->begin(spatial_dimension)[index], true); contact_map[slave_node].projection = Vector(natural_projections->begin(spatial_dimension - 1)[index], true); contact_map[slave_node].connectivity = elem_conn; } } /* -------------------------------------------------------------------------- */ -void ContactDetector::constructGrid(SpatialGrid & grid, BBox & bbox, + void ContactDetector::constructGrid(SpatialGrid & grid, BBox & bbox, const Array & nodes_list) { auto to_grid = [&](UInt node) { Vector pos(spatial_dimension); for (UInt s: arange(spatial_dimension)) { pos(s) = this->positions(node, s); } if (bbox.contains(pos)) { grid.insert(node, pos); } }; std::for_each(nodes_list.begin(), nodes_list.end(), to_grid); } /* -------------------------------------------------------------------------- */ void ContactDetector::constructBoundingBox(BBox & bbox, const Array & nodes_list) { auto to_bbox = [&](UInt node) { Vector pos(spatial_dimension); for (UInt s: arange(spatial_dimension)) { pos(s) = this->positions(node, s); } bbox += pos; }; std::for_each(nodes_list.begin(), nodes_list.end(), to_bbox); AKANTU_DEBUG_INFO("BBox" << bbox); + + std::cerr << bbox << std::endl; auto & lower_bound = bbox.getLowerBounds(); auto & upper_bound = bbox.getUpperBounds(); for (UInt s: arange(spatial_dimension)) { lower_bound(s) -= this->max_bb; upper_bound(s) += this->max_bb; } AKANTU_DEBUG_INFO("BBox" << bbox); } /* -------------------------------------------------------------------------- */ void ContactDetector::computeCellSpacing(Vector & spacing) { for (UInt s: arange(spatial_dimension)) spacing(s) = std::sqrt(2.0) * max_dd; } /* -------------------------------------------------------------------------- */ void ContactDetector::computeOrthogonalProjection(const UInt & node, const Array & elements, Array & normals, Array & gaps, - Array & natural_projections) { + Array & natural_projections, + Array & status) { Vector query(spatial_dimension); for (UInt s: arange(spatial_dimension)) { query(s) = this->positions(node, s); } + std::cerr << "node = " << node << std::endl; for (auto && values : zip( elements, gaps, make_view(normals , spatial_dimension), - make_view(natural_projections, spatial_dimension - 1))) { + make_view(natural_projections, spatial_dimension - 1), + status)) { const auto & element = std::get<0>(values); auto & gap = std::get<1>(values); auto & normal = std::get<2>(values); auto & natural_projection = std::get<3>(values); + auto & to_consider = std::get<4>(values); this->computeNormalOnElement(element, normal); Vector real_projection(spatial_dimension); this->computeProjectionOnElement(element, normal, query, - natural_projection, real_projection); - + natural_projection, real_projection); + Vector distance(spatial_dimension); distance = query - real_projection; gap = Math::norm(spatial_dimension, distance.storage()); - + Vector direction = distance.normalize(); Real cos_angle = direction.dot(normal); - Real tolerance = 1e-8; - + Real tolerance = 1e-8; if (std::abs(cos_angle - 1) <= tolerance && detection_type == _explicit) { gap *= -1; } + + UInt nb_xi_inside = 0; + Real epsilon = 0.05; + for (auto xi : natural_projection) { + if (xi >= -1.0 -epsilon and xi <= 1.0 + epsilon) { + nb_xi_inside++; + } + } + + if (nb_xi_inside == natural_projection.size()) { + to_consider = true; + } + else { + to_consider = false; + } + + if (!to_consider) { + + std::cerr << normal << std::endl; + std::cerr << real_projection << std::endl; + std::cerr << natural_projection << std::endl; + std::cerr << gap << std::endl; + } } - -} +} + /* -------------------------------------------------------------------------- */ void ContactDetector::computeNormalOnElement(const Element & element, Vector & normal) { Matrix vectors(spatial_dimension, spatial_dimension - 1); this->vectorsAlongElement(element, vectors); switch (this->spatial_dimension) { case 2: { Math::normal2(vectors.storage(), normal.storage()); break; } case 3: { Math::normal3(vectors(0).storage(), vectors(1).storage(), normal.storage()); break; } default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); } } } - /* -------------------------------------------------------------------------- */ void ContactDetector::computeProjectionOnElement(const Element & element, const Vector & normal, const Vector & query, Vector & natural_projection, Vector & real_projection) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); Matrix coords(spatial_dimension, nb_nodes_per_element); this->coordinatesOfElement(element, coords); - - Vector point(coords(0)); - + Vector point(coords(0)); Real alpha = (query - point).dot(normal); real_projection = query - alpha * normal; - // use contains function to check whether projection lies inside - // the element, if yes it is a valid projection otherwise no - // still have to think about what to do if normal exists but - // projection doesnot lie inside the element - - bool validity = this->isValidProjection(element, real_projection, natural_projection); - + this->computeNaturalProjection(element, real_projection, natural_projection); } /* -------------------------------------------------------------------------- */ -bool ContactDetector::isValidProjection(const Element & element, Vector & real_projection, - Vector & natural_projection) { +void ContactDetector::computeNaturalProjection(const Element & element, Vector & real_projection, + Vector & natural_projection) { const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(this->positions /*mesh.getNodes()*/, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); #define GET_NATURAL_COORDINATE(type) \ ElementClass::inverseMap(real_projection, nodes_coord, natural_projection) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_COORDINATE); #undef GET_NATURAL_COORDINATE - - /*Vector barycenter(spatial_dimension); - mesh.getBarycenter(element, barycenter); - - Real distance = 0; - switch (this->spatial_dimension) { - case 2: { - distance = Math::distance_2d(real_projection.storage(), barycenter.storage()); - break; - } - case 3: { - distance = Math::distance_3d(real_projection.storage(), barycenter.storage()); - break; - } - default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); } - } - - if (distance <= max_dd) { - return true; - }*/ - - return false; } /* -------------------------------------------------------------------------- */ void ContactDetector::coordinatesOfElement(const Element & el, Matrix & coords) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); Vector connect = mesh.getConnectivity(el.type, _not_ghost) .begin(nb_nodes_per_element)[el.element]; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; for (UInt s: arange(spatial_dimension)) { coords(s, n) = this->positions(node, s); } } } /* -------------------------------------------------------------------------- */ void ContactDetector::vectorsAlongElement(const Element & el, Matrix & vectors) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); Matrix coords(spatial_dimension, nb_nodes_per_element); this->coordinatesOfElement(el, coords); switch (spatial_dimension) { case 2: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); break; } case 3: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); vectors(1) = Vector(coords(2)) - Vector(coords(0)); break; } default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); } } } /* -------------------------------------------------------------------------- */ -void ContactDetector::normalProjection(const Element & el, const Vector & slave_coord, - Vector & natural_coord, Real & tolerance) { +void ContactDetector::normalProjection(const Element & /*el*/, const Vector & /*slave_coord*/, + Vector & /*natural_coord*/, Real & /*tolerance*/) { /*Real fmin; auto update_fmin = [&fmin, &slave_coord, &node_coords, &natural_coord]() { Vector physical_guess_v(physical_guess.storage(), spatial_dimension); // interpolate on natual coordinate and get physical guess // compute gradient or jacobian on natural cooordiante // f = slave_coord - physical_guess; return fmin; }; auto closet_point_error = update_fmin(); while (tolerance < closet_point_error) { // compute gradiend on natural coordinate // compute second variation of shape function at natural coord closet_point_error = update_fmin(); }*/ } } // akantu diff --git a/src/model/contact_mechanics/contact_detector.hh b/src/model/contact_mechanics/contact_detector.hh index 187559b11..a264f11df 100644 --- a/src/model/contact_mechanics/contact_detector.hh +++ b/src/model/contact_mechanics/contact_detector.hh @@ -1,181 +1,188 @@ /** * @file contact_detection.hh * * @author Mohit Pundir * * @date creation: Wed Sep 12 2018 * @date last modification: Tue Oct 23 2018 * * @brief Mother class for all detection algorithms * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_grid_dynamic.hh" #include "aka_bbox.hh" #include "mesh.hh" #include "mesh_io.hh" #include "fe_engine.hh" #include "parsable.hh" #include "element_group.hh" #include "contact_element.hh" #include "element_class.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_DETECTOR_HH__ #define __AKANTU_CONTACT_DETECTOR_HH__ namespace akantu { enum class Surface { master, slave }; /* -------------------------------------------------------------------------- */ class ContactDetector : private Memory, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ public: ContactDetector(Mesh &, const ID & id = "contact_detector", UInt memory_id = 0); ContactDetector(Mesh &, Array & positions, const ID & id = "contact_detector", UInt memory_id = 0); ~ContactDetector() = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ public: /// void search(std::map &); /// computes orthogonal projection on master elements void computeOrthogonalProjection(const UInt & /* slave node */, const Array & /* master elements */, Array & /* normals */, Array & /* gaps */, - Array & /* projections */); + Array & /* projections */, + Array & /* status + */ ); + /// inline std::string getSurfaceId(const std::string name) { if (name == "slave") { return surfaces[Surface::slave]; } }; private: /// reads the input file to get contact detection options void parseSection(); /// performs global spatial search void globalSearch(std::map &); /// performs local search to create contact element /// TODO: templated function typename void localSearch(SpatialGrid &, SpatialGrid &, std::map &); /// constructs a grid containing nodes lying within bounding box /// TODO : templated fucntion to created template Spatial Grid void constructGrid(SpatialGrid &, BBox &, const Array &); /// constructs the bounding box based on nodes list void constructBoundingBox(BBox &, const Array &); /// computes the optimal cell size for grid void computeCellSpacing(Vector &); /// computes the maximum in radius for a given mesh void getMaximalDetectionDistance(); /// extracts the coordinates of an element void coordinatesOfElement(const Element & /* element id */, Matrix & /* coordinates */); /// extracts vectors which forms the plane of element void vectorsAlongElement(const Element & /* element id */, Matrix & /* vectors matrix */); /// computes normal on an element void computeNormalOnElement(const Element & /* element id */, Vector & /* normal vector */); /// computes projection of a query point on an element void computeProjectionOnElement(const Element & /* element */, const Vector & /* normal */, const Vector & /* query */, Vector & /* projection */, Vector & /* real_projection */); - /// checks for the validity of a projection - bool isValidProjection(const Element & /* element */, - Vector & /* real projection */, - Vector & /* natural projection */); + + /// computes natural projection of a real projection + void computeNaturalProjection(const Element & /* element */, + Vector & /* real projection */, + Vector & /* natural projection */); /// compute normal projection of slav coord on a given element void normalProjection(const Element & el, const Vector & slave_coord, Vector & natural_coord, Real & tolerance); + /// compute the active zones for implicit detection type, uses + /// displacement increment to compute the actual gap value for slave + /// nodes which will penetrate the master surface /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// maximal detection distance for grid spacing Real max_dd; /// maximal bounding box extension Real max_bb; /// Mesh Mesh & mesh; /// dimension of the model UInt spatial_dimension{0}; /// map to contain ids for surfaces std::map surfaces; /// contains the updated positions of the nodes Array & positions; /// type of detection extrinisic/intrinsic ContactDetectorType detection_type; }; } // namespace akantu #endif /* __AKANTU_CONTACT_DETECTOR_HH__ */ diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index a1f87b03d..bcb54e356 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,588 +1,659 @@ /** * @file coontact_mechanics_model.cc * * @author Mohit Pundir * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Contact mechanics model * * @section LICENSE * * Copyright (©) 2010-2018 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 "contact_mechanics_model.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id), current_positions(mesh.getNodes()) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->initDOFManager(); this->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, id + ":contact_detector"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, Array & positions, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id), current_positions(positions) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->initDOFManager(); this->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, positions, id + ":contact_detector"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::~ContactMechanicsModel() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // initalize the resolutions if (this->parser.getLastParsedFile() != "") { this->instantiateResolutions(); } this->initResolutions(); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::instantiateResolutions() { ParserSection model_section; bool is_empty; std::tie(model_section, is_empty) = this->getParserSection(); if (not is_empty) { auto model_resolutions = model_section.getSubSections(ParserType::_contact_resolution); for (const auto & section : model_resolutions) { this->registerNewResolution(section); } } auto sub_sections = this->parser.getSubSections(ParserType::_contact_resolution); for (const auto & section : sub_sections) { this->registerNewResolution(section); } if (resolutions.empty()) AKANTU_EXCEPTION("No contact resolutions where instantiated for the model" << getID()); are_resolutions_instantiated = true; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution(const ParserSection & section) { std::string res_name; std::string res_type = section.getName(); std::string opt_param = section.getOption(); try { std::string tmp = section.getParameter("name"); res_name = tmp; /** this can seem weird, but there is an ambiguous operator * overload that i couldn't solve. @todo remove the * weirdness of this code */ } catch (debug::Exception &) { AKANTU_ERROR("A contact resolution of type \'" << res_type << "\' in the input file has been defined without a name!"); } Resolution & res = this->registerNewResolution(res_name, res_type, opt_param); res.parseSection(section); return res; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution(const ID & res_name, const ID & res_type, const ID & opt_param) { AKANTU_DEBUG_ASSERT(resolutions_names_to_id.find(res_name) == resolutions_names_to_id.end(), "A resolution with this name '" << res_name << "' has already been registered. " << "Please use unique names for resolutions"); UInt res_count = resolutions.size(); resolutions_names_to_id[res_name] = res_count; std::stringstream sstr_res; sstr_res << this->id << ":" << res_count << ":" << res_type; ID res_id = sstr_res.str(); std::unique_ptr resolution = ResolutionFactory::getInstance().allocate( res_type, spatial_dimension, opt_param, *this, res_id); resolutions.push_back(std::move(resolution)); return *(resolutions.back()); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initResolutions() { AKANTU_DEBUG_ASSERT(resolutions.size() != 0, "No resolutions to initialize !"); if (!are_resolutions_instantiated) instantiateResolutions(); // \TODO check if each resolution needs a initResolution() method } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->contact_force, spatial_dimension, "contact_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->areas, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->contact_force, _dst_nodal); //dof_manager.registerDOFs("displacement", *this->blocked_dofs); } } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", _tsst_dynamic); } case _implicit_contact: { return std::make_tuple("implicit_contact", _tsst_static); } default: return std::make_tuple("unkown", _tsst_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type ) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the contact forces"); // assemble the forces due to local stresses auto assemble = [&] (auto && ghost_type) { for (auto & resolution : resolutions) { resolution->assembleInternalForces(ghost_type); } }; AKANTU_DEBUG_INFO("Assemble residual for local elements"); assemble(_not_ghost); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); assemble(_ghost); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->contact_force, 1); AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { this->detector->search(this->contact_map); auto & blocked_dof = const_cast &>(this->getBlockedDOFs()); auto & gap = const_cast &>(this->getGaps()); for(auto & entry : contact_map) { const auto & element = entry.second; const auto & connectivity = element.connectivity; for (UInt i = 0; i < connectivity.size(); ++i) { UInt n = connectivity(i); blocked_dof[n] = 1.0; gap[n] = element.gap; } } this->areas->clear(); this->external_force->clear(); this->applyBC(BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), this->detector->getSurfaceId("slave")); auto ext_force_it = external_force->begin(Model::spatial_dimension); auto areas_it = areas->begin(); UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++areas_it) { const auto & ext_force = *ext_force_it; auto & area = *areas_it; for (UInt i = 0; i < Model::spatial_dimension; ++i) { area += pow(ext_force(i), 2); } area = sqrt(area); } } + + +/* -------------------------------------------------------------------------- */ +void ContactMechanicsModel::search(Array & increment) { + + this->detector->search(this->contact_map); + + auto & blocked_dof = + const_cast &>(this->getBlockedDOFs()); + + auto & gap = + const_cast &>(this->getGaps()); + + for (auto & entry: contact_map) { + auto & element = entry.second; + const auto & connectivity = element.connectivity; + + auto master_node = connectivity[1]; + Vector u(spatial_dimension); + for (UInt s : arange(spatial_dimension)) { + u(s) = increment(master_node, s); + } + + u *= -1.0; + std::cerr << u << std::endl; + + const auto & normal = element.normal; + std::cerr << normal << std::endl; + std::cerr << element.gap << std::endl; + Real uv = Math::vectorDot(u.storage(), normal.storage(), spatial_dimension); + + if (uv + element.gap <= 0) { + element.gap = abs(uv + element.gap); + } + else { + element.gap = 0.0; + } + } + + + for(auto & entry : contact_map) { + const auto & element = entry.second; + const auto & connectivity = element.connectivity; + for (UInt i = 0; i < connectivity.size(); ++i) { + UInt n = connectivity(i); + blocked_dof[n] = 1.0; + gap[n] = element.gap; + } + } + + this->areas->clear(); + this->external_force->clear(); + + this->applyBC(BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), + this->detector->getSurfaceId("slave")); + + auto ext_force_it = external_force->begin(Model::spatial_dimension); + auto areas_it = areas->begin(); + UInt nb_nodes = this->mesh.getNbNodes(); + + for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++areas_it) { + const auto & ext_force = *ext_force_it; + auto & area = *areas_it; + + for (UInt i = 0; i < Model::spatial_dimension; ++i) { + area += pow(ext_force(i), 2); + } + + area = sqrt(area); + } +} + /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() { } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep() { } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Contact Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + resolutions [" << std::endl; for (auto & resolution : resolutions) { resolution->printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ MatrixType ContactMechanicsModel::getMatrixType(const ID & matrix_id) { // \TODO correct it for contact mechanics model, only one type of matrix if (matrix_id == "K") return _symmetric; return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().clearMatrix("K"); for (auto & resolution : resolutions) { resolution->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * ContactMechanicsModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { //std::map *> uint_nodal_fields; //uint_nodal_fields["blocked_dofs"] = blocked_dofs; dumper::Field * field = nullptr; //field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * ContactMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["contact_force"] = this->contact_force; real_nodal_fields["external_force"] = this->external_force; real_nodal_fields["blocked_dofs"] = this->blocked_dofs; real_nodal_fields["gaps"] = this->gaps; real_nodal_fields["areas"] = this->areas; dumper::Field * field = nullptr; if (padding_flag) field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); else field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } #else /* -------------------------------------------------------------------------- */ dumper::Field * ContactMechanicsModel::createNodalFieldReal( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData(const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/contact_mechanics/contact_mechanics_model.hh b/src/model/contact_mechanics/contact_mechanics_model.hh index b192efd7a..4213ba4f7 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,317 +1,319 @@ /** * @file contact_mechanics_model.hh * * @author Mohit Pundir * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Contact Mechanics * * @section LICENSE * * Copyright (©) 2010-2018 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 "boundary_condition.hh" #include "data_accessor.hh" #include "model.hh" #include "fe_engine.hh" #include "contact_detector.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_MECHANICS_MODEL_HH__ #define __AKANTU_CONTACT_MECHANICS_MODEL_HH__ namespace akantu { class Resolution; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class ContactMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: ContactMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_contact_mechanics_model); ContactMechanicsModel(Mesh & mesh, Array & positions, UInt spatial_dimension = _all_dimensions, const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_contact_mechanics_model); ~ContactMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// allocate all vectors void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize all internal arrays for resolutions void initResolutions(); /// initialize the modelType void initModel() override; /// call back for the solver, computes the force residual void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Contact Detection */ /* ------------------------------------------------------------------------ */ public: void search(); + + void search(Array &); template void computeNodalAreas(const FunctorType &, const std::string &, GhostType); /* ------------------------------------------------------------------------ */ /* Contact Resolution */ /* ------------------------------------------------------------------------ */ public: /// register an empty contact resolution of a given type Resolution & registerNewResolution(const ID & res_name, const ID & res_type, const ID & opt_param); protected: /// register a resolution in the dynamic database Resolution & registerNewResolution(const ParserSection & res_section); /// read the resolution files to instantiate all the resolutions void instantiateResolutions(); /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: /// contact detection class friend class ContactDetector; /// contact resolution class friend class Resolution; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the ContactMechanics::contact_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *contact_force, Array &); /// get the SolidMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the SolidMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the ContactMechanics::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the ContactMechanics::gaps vector (contact gaps) AKANTU_GET_MACRO(Gaps, *gaps, Array &); /// get the ContactMechanics::areas vector (contact areas) AKANTU_GET_MACRO(ContactArea, *areas, Array &); /// get the ContactMechanics::contact_force vector (internal forces) AKANTU_GET_MACRO(CurrentPositions, current_positions, Array &); /// get the contat map inline std::map & getContactMap() { return contact_map; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// tells if the resolutions are instantiated bool are_resolutions_instantiated; /// displacements array Array * displacement; /// increment of displacement Array * displacement_increment{nullptr}; /// contact forces array Array * contact_force{nullptr}; /// external forces array Array * external_force{nullptr}; /// boundary vector Array * blocked_dofs{nullptr}; /// gaps vector Array * gaps{nullptr}; /// nodal areas vector Array * areas{nullptr}; /// array of current position used during update residual Array & current_positions; /// contact detection std::unique_ptr detector; /// list of contact resolutions std::vector> resolutions; /// mapping between resolution name and resolution internal id std::map resolutions_names_to_id; /// mapping between slave node its respective contact element std::map contact_map; }; } // namespace akantu /* ------------------------------------------------------------------------ */ /* inline functions */ /* ------------------------------------------------------------------------ */ #include "resolution.hh" #include "parser.hh" /* ------------------------------------------------------------------------ */ #endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc index 81f90cbc0..0d0bdb889 100644 --- a/src/model/contact_mechanics/resolution.cc +++ b/src/model/contact_mechanics/resolution.cc @@ -1,363 +1,384 @@ /** * @file resolution.cc * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Implementation of common part of the contact resolution class * * @section LICENSE * * Copyright (©) 2010-2018 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 "resolution.hh" #include "contact_mechanics_model.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Resolution::Resolution(ContactMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_contact_resolution, id), fem(model.getFEEngine()), name(""), model(model), spatial_dimension(model.getMesh().getSpatialDimension()){ AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Resolution::~Resolution() = default; /* -------------------------------------------------------------------------- */ void Resolution::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("mu", mu, Real(0.), _pat_parsable | _pat_modifiable, "Friciton Coefficient"); } /* -------------------------------------------------------------------------- */ void Resolution::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "Contact Resolution " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Resolution::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & internal_force = const_cast &>(model.getInternalForce()); auto & contact_area = const_cast &>(model.getContactArea()); auto & contact_map = model.getContactMap(); const auto slave_nodes = model.getMesh().getElementGroup(name).getNodes(); for (auto & slave: slave_nodes) { if (contact_map.find(slave) == contact_map.end()) continue; auto & master = contact_map[slave].master; auto & gap = contact_map[slave].gap; auto & projection = contact_map[slave].projection; auto & normal = contact_map[slave].normal; const auto & connectivity = contact_map[slave].connectivity; const ElementType & type = master.type; UInt nb_nodes_master = Mesh::getNbNodesPerElement(master.type); Vector shapes(nb_nodes_master); Matrix shapes_derivatives(spatial_dimension - 1, nb_nodes_master); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shapes_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Vector elem_force(connectivity.size() * spatial_dimension); Matrix tangents(spatial_dimension - 1, spatial_dimension); Matrix global_coords(nb_nodes_master, spatial_dimension); computeCoordinates(master, global_coords); computeTangents(shapes_derivatives, global_coords, tangents); - + Matrix surface_matrix(spatial_dimension - 1, spatial_dimension - 1); computeSurfaceMatrix(tangents, surface_matrix); Vector n(connectivity.size() * spatial_dimension); computeN(n, shapes, normal); computeNormalForce(elem_force, n, gap); Array t_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); Array n_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); Array d_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); computeTalpha(t_alpha, shapes, tangents); computeNalpha(n_alpha, shapes_derivatives, normal); computeDalpha(d_alpha, n_alpha, t_alpha, surface_matrix, gap); //computeFrictionForce(elem_force, d_alpha, gap); UInt nb_degree_of_freedom = internal_force.getNbComponent(); for (UInt i = 0; i < connectivity.size(); ++i) { UInt n = connectivity[i]; for (UInt j = 0; j < nb_degree_of_freedom; ++j) { UInt offset_node = n * nb_degree_of_freedom + j; internal_force[offset_node] += elem_force[i*nb_degree_of_freedom + j]; internal_force[offset_node] *= contact_area[n]; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & contact_stiffness = const_cast(model.getDOFManager().getMatrix("K")); const auto slave_nodes = model.getMesh().getElementGroup(name).getNodes(); + auto & contact_area = + const_cast &>(model.getContactArea()); + auto & contact_map = model.getContactMap(); for (auto & slave: slave_nodes) { if (contact_map.find(slave) == contact_map.end()) { continue; } auto & master = contact_map[slave].master; auto & gap = contact_map[slave].gap; auto & projection = contact_map[slave].projection; auto & normal = contact_map[slave].normal; const auto & connectivity = contact_map[slave].connectivity; const ElementType & type = master.type; UInt nb_nodes_master = Mesh::getNbNodesPerElement(master.type); Vector shapes(nb_nodes_master); Matrix shapes_derivatives(spatial_dimension - 1, nb_nodes_master); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shapes_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Matrix elementary_stiffness(connectivity.size() * spatial_dimension, connectivity.size() * spatial_dimension); Matrix tangents(spatial_dimension - 1, spatial_dimension); Matrix global_coords(nb_nodes_master, spatial_dimension); computeCoordinates(master, global_coords); computeTangents(shapes_derivatives, global_coords, tangents); + std::cerr << "slave" << slave << std::endl; + std::cerr << "normal = " << normal << std::endl; + std::cerr << "projection = " << projection << std::endl; + std::cerr << "coords = " << global_coords << std::endl; + std::cerr << "shapes = "<< shapes << std::endl; + std::cerr << "derivatives = " << shapes_derivatives << std::endl; + std::cerr << "tangents = " << tangents << std::endl; + Matrix surface_matrix(spatial_dimension - 1, spatial_dimension - 1); computeSurfaceMatrix(tangents, surface_matrix); Vector n(connectivity.size() * spatial_dimension); Array t_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); Array n_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); Array d_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); - + computeN( n, shapes, normal); computeTalpha( t_alpha, shapes, tangents); computeNalpha( n_alpha, shapes_derivatives, normal); computeDalpha( d_alpha, n_alpha, t_alpha, surface_matrix, gap); - + Matrix kc(connectivity.size() * spatial_dimension, connectivity.size() * spatial_dimension); computeTangentModuli(kc, n, n_alpha, d_alpha, surface_matrix, gap); - + std::vector equations; UInt nb_degree_of_freedom = model.getSpatialDimension(); + + std::vector areas; for (UInt i : arange(connectivity.size())) { UInt n = connectivity[i]; - for (UInt j : arange(nb_degree_of_freedom)) + for (UInt j : arange(nb_degree_of_freedom)) { equations.push_back(n * nb_degree_of_freedom + j); + areas.push_back(contact_area[n]); + } + } + + /*for (auto c : connectivity) { + std::cerr << c << " "; } + std::cerr << "\n"; + + for (auto e : equations) { + std::cerr << e << " "; + } + std::cerr << "\n";*/ + for (UInt i : arange(kc.rows())) { UInt row = equations[i]; for (UInt j : arange(kc.cols())) { UInt col = equations[j]; + //std::cerr << row << "," << col << "=" << kc(i, j) < & shapes_derivatives, Matrix & global_coords, Matrix & tangents) { - /*UInt index = 0; - for (auto && values : zip(make_view(tangents, spatial_dimension))) { - auto & tangent = std::get<0>(values); - for (UInt n : arange(global_coords.getNbComponent())) { - tangent += shapes_derivatives(n, index) * global_coords(n); - } - ++index; - }*/ - tangents.mul(shapes_derivatives, global_coords); + tangents.mul(shapes_derivatives, global_coords); } /* -------------------------------------------------------------------------- */ void Resolution::computeSurfaceMatrix(Matrix & tangents, Matrix & surface_matrix) { - /*Matrix A(surface_matrix); - for (UInt i : arange(spatial_dimension - 1)) { - for (UInt j : arange(spatial_dimension -1 )) { - A(i, j) = tangents(i) * tangents(j); - } - }*/ - surface_matrix.mul(tangents, tangents); + //std::cerr << "surface = " << surface_matrix << std::endl; surface_matrix = surface_matrix.inverse(); - //A.inverse(surface_matrix); } /* -------------------------------------------------------------------------- */ void Resolution::computeN(Vector & n, Vector & shapes, Vector & normal) { UInt dim = normal.size(); for (UInt i = 0; i < dim; ++i) { n[i] = normal[i]; for (UInt j = 0; j < shapes.size(); ++j) { n[(1 + j) * dim + i] = -normal[i] * shapes[j]; } } } /* -------------------------------------------------------------------------- */ void Resolution::computeTalpha(Array & t_alpha, Vector & shapes, Matrix & tangents) { - + t_alpha.clear(); for (auto && values: zip(tangents.transpose(), make_view(t_alpha, t_alpha.size()))) { auto & tangent = std::get<0>(values); auto & t_s = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { t_s[i] = -tangent(i); for (UInt j : arange(shapes.size())) { t_s[(1 + j)*spatial_dimension + i] = -shapes[j] * tangent(i); } } + //std::cerr << "t_s = " << t_s << std::endl; } } /* -------------------------------------------------------------------------- */ void Resolution::computeNalpha(Array & n_alpha, Matrix & shapes_derivatives, Vector & normal) { - + n_alpha.clear(); for (auto && values: zip(shapes_derivatives.transpose(), make_view(n_alpha, n_alpha.size()))) { auto & dnds = std::get<0>(values); auto & n_s = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { n_s[i] = 0; - for (UInt j : arange(shapes_derivatives.size())) { - n_s[(1 + j)*spatial_dimension + i] = -shapes_derivatives[j]*normal[i]; + for (UInt j : arange(dnds.size())) { + n_s[(1 + j)*spatial_dimension + i] = -dnds(j)*normal[i]; } } + //std::cerr << "n_s = " << n_s << std::endl; } } /* -------------------------------------------------------------------------- */ void Resolution::computeDalpha(Array & d_alpha, Array & n_alpha, Array & t_alpha, Matrix & surface_matrix, Real & gap) { + d_alpha.clear(); for (auto && entry : zip(surface_matrix.transpose(), make_view(d_alpha, d_alpha.size()))) { auto & a_s = std::get<0>(entry); auto & d_s = std::get<1>(entry); for (auto && values : zip(arange(t_alpha.size()), make_view(t_alpha, t_alpha.size()), make_view(n_alpha, n_alpha.size()))) { auto & index = std::get<0>(values); auto & t_s = std::get<1>(values); auto & n_s = std::get<2>(values); + //std::cerr << "d_s = " << d_s << std::endl; d_s += (t_s + gap * n_s); + //std::cerr << "d_s = " << d_s << std::endl; d_s *= a_s(index); } + //std::cerr << "d_s" << d_s << std::endl; } } /* -------------------------------------------------------------------------- */ void Resolution::computeCoordinates(const Element & el, Matrix & coords) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); Vector connect = model.getMesh().getConnectivity(el.type, _not_ghost) .begin(nb_nodes_per_element)[el.element]; // change this to current position auto & positions = model.getMesh().getNodes(); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; for (UInt s: arange(spatial_dimension)) { coords(n, s) = positions(node, s); } } } } // akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc index 3a280bc3a..cfd6452e9 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc @@ -1,160 +1,165 @@ /** * @file resolution_penalty.cc * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Specialization of the resolution class for the penalty method * * @section LICENSE * * Copyright (©) 2010-2018 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 "resolution_penalty.hh" namespace akantu { /* -------------------------------------------------------------------------- */ ResolutionPenalty::ResolutionPenalty(ContactMechanicsModel & model, const ID & id) : Resolution(model, id) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::initialize() { this->registerParam("epsilon", epsilon, Real(0.), _pat_parsable | _pat_modifiable, "Penalty parameter"); } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeNormalForce(Vector & force, Vector & n, Real & gap) { + + force.clear(); Real tn = gap * epsilon; tn = macaulay(tn); for (UInt i : arange(force.size())) { force[i] += tn * n[i]; } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeFrictionForce(Vector & force, Array & d_alpha, Real & gap) { Vector tractions(d_alpha.getNbComponent()); computeFrictionalTraction(tractions); for (auto && values: zip(tractions, make_view(d_alpha, d_alpha.size()))) { auto & t_s = std::get<0>(values); auto & d_s = std::get<1>(values); force += d_s * t_s; } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeFrictionalTraction(Vector & tractions) { } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTangentModuli(Matrix & kc, Vector & n, Array & n_alpha, Array & d_alpha, Matrix & surface_matrix, Real & gap) { computeNormalStiffness(kc, n, n_alpha, d_alpha, surface_matrix, gap); computeFrictionalStiffness(n, n_alpha, d_alpha, gap); } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeNormalStiffness(Matrix & kc, Vector & n, +void ResolutionPenalty::computeNormalStiffness(Matrix & ke, Vector & n, Array & n_alpha, Array & d_alpha, Matrix & surface_matrix, Real & gap) { Real tn = gap * epsilon; tn = macaulay(tn); - Matrix ke(n_alpha.size(), n_alpha.size()); + //Matrix ke(n_alpha.size(), n_alpha.size()); Matrix n_mat(n.storage(), n.size(), 1); ke.mul(n_mat, n_mat); ke *= epsilon * heaviside(gap); + + std::cerr << " ke = " << ke << std::endl; - for (auto && entry1 : + /*for (auto && entry1 : enumerate(make_view(n_alpha, n_alpha.size()))) { auto & i = std::get<0>(entry1); auto & ni = std::get<1>(entry1); Matrix ni_mat(ni.storage(), ni.size(), 1); for (auto && entry2: enumerate(make_view(n_alpha, n_alpha.size()))) { auto & j = std::get<0>(entry2); auto & nj = std::get<1>(entry2); Matrix nj_mat(nj.storage(), nj.size(), 1); Matrix tmp(nj.size(), nj.size()); tmp.mul(ni_mat, nj_mat); - tmp *= surface_matrix[i, j]; + tmp *= surface_matrix[i, j] * tn; ke += tmp; } - } + }*/ for (auto && values: zip(make_view(n_alpha, n_alpha.size()), make_view(d_alpha, d_alpha.size()))) { auto & n_s = std::get<0>(values); auto & d_s = std::get<1>(values); Matrix ns_mat(n_s.storage(), n_s.size(), 1); Matrix ds_mat(d_s.storage(), d_s.size(), 1); Matrix tmp1(n_s.size(), n_s.size()); tmp1.mul(ns_mat, ds_mat); Matrix tmp2(n_s.size(), n_s.size()); tmp1.mul(ds_mat, ns_mat); - ke -= tmp1 + tmp2; + ke -= (tmp1 + tmp2) * tn; } - + + std::cerr << " ke = " << ke << std::endl; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeFrictionalStiffness(Vector & n, Array & n_alpha, Array & d_alpha, Real & gap) { } INSTANTIATE_RESOLUTION(penalty, ResolutionPenalty); } // akantu diff --git a/src/model/integration_scheme/pseudo_time.cc b/src/model/integration_scheme/pseudo_time.cc index 1ca15129b..93f9a40fb 100644 --- a/src/model/integration_scheme/pseudo_time.cc +++ b/src/model/integration_scheme/pseudo_time.cc @@ -1,82 +1,82 @@ /** * @file pseudo_time.cc * * @author Nicolas Richart * * @date creation: Fri Feb 19 2016 * @date last modification: Wed Jan 31 2018 * * @brief Implementation of a really simple integration scheme * * @section LICENSE * * Copyright (©) 2016-2018 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 "pseudo_time.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ PseudoTime::PseudoTime(DOFManager & dof_manager, const ID & dof_id) : IntegrationScheme(dof_manager, dof_id, 0), k_release(0) {} /* -------------------------------------------------------------------------- */ std::vector PseudoTime::getNeededMatrixList() { return {"K"}; } /* -------------------------------------------------------------------------- */ void PseudoTime::predictor(Real) {} /* -------------------------------------------------------------------------- */ void PseudoTime::corrector(const SolutionType &, Real) { auto & us = this->dof_manager.getDOFs(this->dof_id); const auto & deltas = this->dof_manager.getSolution(this->dof_id); const auto & blocked_dofs = this->dof_manager.getBlockedDOFs(this->dof_id); for (auto && tuple : zip(make_view(us), deltas, make_view(blocked_dofs))) { auto & u = std::get<0>(tuple); const auto & delta = std::get<1>(tuple); const auto & bld = std::get<2>(tuple); if (not bld) u += delta; } } /* -------------------------------------------------------------------------- */ void PseudoTime::assembleJacobian(const SolutionType &, Real) { SparseMatrix & J = this->dof_manager.getMatrix("J"); const SparseMatrix & K = this->dof_manager.getMatrix("K"); if (K.getRelease() == k_release) return; J.clear(); J.add(K); - + J.saveMatrix("J.mtx"); k_release = K.getRelease(); } /* -------------------------------------------------------------------------- */ void PseudoTime::assembleResidual(bool) {} /* -------------------------------------------------------------------------- */ } // akantu diff --git a/src/model/model_solver.cc b/src/model/model_solver.cc index 1891e7fd3..5ef247a5a 100644 --- a/src/model/model_solver.cc +++ b/src/model/model_solver.cc @@ -1,365 +1,371 @@ /** * @file model_solver.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of ModelSolver * * @section LICENSE * * Copyright (©) 2015-2018 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 "model_solver.hh" #include "dof_manager.hh" #include "dof_manager_default.hh" #include "mesh.hh" #include "non_linear_solver.hh" #include "time_step_solver.hh" #if defined(AKANTU_USE_PETSC) #include "dof_manager_petsc.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template static T getOptionToType(const std::string & opt_str) { std::stringstream sstr(opt_str); T opt; sstr >> opt; return opt; } /* -------------------------------------------------------------------------- */ ModelSolver::ModelSolver(Mesh & mesh, const ModelType & type, const ID & id, UInt memory_id) : Parsable(ParserType::_model, id), SolverCallback(), model_type(type), parent_id(id), parent_memory_id(memory_id), mesh(mesh), dof_manager(nullptr), default_solver_id("") {} /* -------------------------------------------------------------------------- */ ModelSolver::~ModelSolver() = default; /* -------------------------------------------------------------------------- */ std::tuple ModelSolver::getParserSection() { auto sub_sections = getStaticParser().getSubSections(ParserType::_model); auto it = std::find_if( sub_sections.begin(), sub_sections.end(), [&](auto && section) { ModelType type = getOptionToType(section.getName()); // default id should be the model type if not defined std::string name = section.getParameter("name", this->parent_id); return type == model_type and name == this->parent_id; }); if (it == sub_sections.end()) return std::make_tuple(ParserSection(), true); return std::make_tuple(*it, false); } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager() { // default without external solver activated at compilation same as mumps that // is the historical solver but with only the lumped solver ID solver_type = "default"; #if defined(AKANTU_USE_MUMPS) solver_type = "default"; #elif defined(AKANTU_USE_PETSC) solver_type = "petsc"; #endif ParserSection section; bool is_empty; std::tie(section, is_empty) = this->getParserSection(); if (not is_empty) { solver_type = section.getOption(solver_type); this->initDOFManager(section, solver_type); } else { this->initDOFManager(solver_type); } } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager(const ID & solver_type) { try { this->dof_manager = DOFManagerFactory::getInstance().allocate( solver_type, mesh, this->parent_id + ":dof_manager" + solver_type, this->parent_memory_id); } catch (...) { AKANTU_EXCEPTION( "To use the solver " << solver_type << " you will have to code it. This is an unknown solver type."); } this->setDOFManager(*this->dof_manager); } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager(const ParserSection & section, const ID & solver_type) { this->initDOFManager(solver_type); auto sub_sections = section.getSubSections(ParserType::_time_step_solver); // parsing the time step solvers for (auto && section : sub_sections) { ID type = section.getName(); ID solver_id = section.getParameter("name", type); auto tss_type = getOptionToType(type); auto tss_options = this->getDefaultSolverOptions(tss_type); auto sub_solvers_sect = section.getSubSections(ParserType::_non_linear_solver); auto nb_non_linear_solver_section = section.getNbSubSections(ParserType::_non_linear_solver); auto nls_type = tss_options.non_linear_solver_type; if (nb_non_linear_solver_section == 1) { auto && nls_section = *(sub_solvers_sect.first); nls_type = getOptionToType(nls_section.getName()); } else if (nb_non_linear_solver_section > 0) { AKANTU_EXCEPTION("More than one non linear solver are provided for the " "time step solver " << solver_id); } this->getNewSolver(solver_id, tss_type, nls_type); if (nb_non_linear_solver_section == 1) { const auto & nls_section = *(sub_solvers_sect.first); this->dof_manager->getNonLinearSolver(solver_id).parseSection( nls_section); } auto sub_integrator_sections = section.getSubSections(ParserType::_integration_scheme); for (auto && is_section : sub_integrator_sections) { const auto & dof_type_str = is_section.getName(); ID dof_id; try { ID tmp = is_section.getParameter("name"); dof_id = tmp; } catch (...) { AKANTU_EXCEPTION("No degree of freedom name specified for the " "integration scheme of type " << dof_type_str); } auto it_type = getOptionToType(dof_type_str); IntegrationScheme::SolutionType s_type = is_section.getParameter( "solution_type", tss_options.solution_type[dof_id]); this->setIntegrationScheme(solver_id, dof_id, it_type, s_type); } for (auto & is_type : tss_options.integration_scheme_type) { if (!this->hasIntegrationScheme(solver_id, is_type.first)) { this->setIntegrationScheme(solver_id, is_type.first, is_type.second, tss_options.solution_type[is_type.first]); } } } if (section.hasParameter("default_solver")) { ID default_solver = section.getParameter("default_solver"); if (this->hasSolver(default_solver)) { this->setDefaultSolver(default_solver); } else AKANTU_EXCEPTION( "The solver \"" << default_solver << "\" was not created, it cannot be set as default solver"); } } /* -------------------------------------------------------------------------- */ TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) { ID tmp_solver_id = solver_id; if (tmp_solver_id == "") tmp_solver_id = this->default_solver_id; TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id); return tss; } /* -------------------------------------------------------------------------- */ const TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) const { ID tmp_solver_id = solver_id; if (solver_id == "") tmp_solver_id = this->default_solver_id; const TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id); return tss; } /* -------------------------------------------------------------------------- */ TimeStepSolver & ModelSolver::getTimeStepSolver(const ID & solver_id) { return this->getSolver(solver_id); } /* -------------------------------------------------------------------------- */ const TimeStepSolver & ModelSolver::getTimeStepSolver(const ID & solver_id) const { return this->getSolver(solver_id); } /* -------------------------------------------------------------------------- */ NonLinearSolver & ModelSolver::getNonLinearSolver(const ID & solver_id) { return this->getSolver(solver_id).getNonLinearSolver(); } /* -------------------------------------------------------------------------- */ const NonLinearSolver & ModelSolver::getNonLinearSolver(const ID & solver_id) const { return this->getSolver(solver_id).getNonLinearSolver(); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasSolver(const ID & solver_id) const { ID tmp_solver_id = solver_id; if (solver_id == "") tmp_solver_id = this->default_solver_id; if (not this->dof_manager) { AKANTU_EXCEPTION("No DOF manager was initialized"); } return this->dof_manager->hasTimeStepSolver(tmp_solver_id); } /* -------------------------------------------------------------------------- */ void ModelSolver::setDefaultSolver(const ID & solver_id) { AKANTU_DEBUG_ASSERT( this->hasSolver(solver_id), "Cannot set the default solver to a solver that does not exists"); this->default_solver_id = solver_id; } /* -------------------------------------------------------------------------- */ void ModelSolver::solveStep(const ID & solver_id) { + solveStep(*this, solver_id); +} + +/* -------------------------------------------------------------------------- */ +void ModelSolver::solveStep(SolverCallback & callback, + const ID & solver_id) { AKANTU_DEBUG_IN(); TimeStepSolver & tss = this->getSolver(solver_id); // make one non linear solve - tss.solveStep(*this); + tss.solveStep(callback); AKANTU_DEBUG_OUT(); } - + /* -------------------------------------------------------------------------- */ void ModelSolver::getNewSolver(const ID & solver_id, TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { if (this->default_solver_id == "") { this->default_solver_id = solver_id; } if (non_linear_solver_type == _nls_auto) { switch (time_step_solver_type) { case _tsst_dynamic: case _tsst_static: non_linear_solver_type = _nls_newton_raphson; break; case _tsst_dynamic_lumped: non_linear_solver_type = _nls_lumped; break; case _tsst_not_defined: AKANTU_EXCEPTION(time_step_solver_type << " is not a valid time step solver type"); break; } } this->initSolver(time_step_solver_type, non_linear_solver_type); NonLinearSolver & nls = this->dof_manager->getNewNonLinearSolver( solver_id, non_linear_solver_type); this->dof_manager->getNewTimeStepSolver(solver_id, time_step_solver_type, nls); } /* -------------------------------------------------------------------------- */ Real ModelSolver::getTimeStep(const ID & solver_id) const { const TimeStepSolver & tss = this->getSolver(solver_id); return tss.getTimeStep(); } /* -------------------------------------------------------------------------- */ void ModelSolver::setTimeStep(Real time_step, const ID & solver_id) { TimeStepSolver & tss = this->getSolver(solver_id); return tss.setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void ModelSolver::setIntegrationScheme( const ID & solver_id, const ID & dof_id, const IntegrationSchemeType & integration_scheme_type, IntegrationScheme::SolutionType solution_type) { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); tss.setIntegrationScheme(dof_id, integration_scheme_type, solution_type); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasDefaultSolver() const { return (this->default_solver_id != ""); } /* -------------------------------------------------------------------------- */ bool ModelSolver::hasIntegrationScheme(const ID & solver_id, const ID & dof_id) const { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); return tss.hasIntegrationScheme(dof_id); } /* -------------------------------------------------------------------------- */ void ModelSolver::predictor() {} /* -------------------------------------------------------------------------- */ void ModelSolver::corrector() {} /* -------------------------------------------------------------------------- */ TimeStepSolverType ModelSolver::getDefaultSolverType() const { return _tsst_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions ModelSolver::getDefaultSolverOptions(__attribute__((unused)) const TimeStepSolverType & type) const { ModelSolverOptions options; options.non_linear_solver_type = _nls_auto; return options; } } // namespace akantu diff --git a/src/model/model_solver.hh b/src/model/model_solver.hh index 7b2c61c21..ab722d33a 100644 --- a/src/model/model_solver.hh +++ b/src/model/model_solver.hh @@ -1,194 +1,201 @@ /** * @file model_solver.hh * * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Wed Jan 31 2018 * * @brief Class regrouping the common solve interface to the different models * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "integration_scheme.hh" #include "parsable.hh" #include "solver_callback.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MODEL_SOLVER_HH__ #define __AKANTU_MODEL_SOLVER_HH__ namespace akantu { class Mesh; class DOFManager; class TimeStepSolver; class NonLinearSolver; struct ModelSolverOptions; } namespace akantu { class ModelSolver : public Parsable, public SolverCallback, public SynchronizerRegistry { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ModelSolver(Mesh & mesh, const ModelType & type, const ID & id, UInt memory_id); ~ModelSolver() override; /// initialize the dof manager based on solver type passed in the input file void initDOFManager(); /// initialize the dof manager based on the used chosen solver type void initDOFManager(const ID & solver_type); protected: /// initialize the dof manager based on the used chosen solver type void initDOFManager(const ParserSection & section, const ID & solver_type); /// Callback for the model to instantiate the matricees when needed virtual void initSolver(TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType /*non_linear_solver_type*/) {} /// get the section in the input file (if it exsits) corresponding to this /// model std::tuple getParserSection(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// solve a step using a given pre instantiated time step solver and - /// nondynamic linear solver + /// non linear solver virtual void solveStep(const ID & solver_id = ""); + /// solve a step using a given pre instantiated time step solver and + /// non linear solver with a user defined callback instead of the + /// model itself /!\ This can mess up everything + virtual void solveStep(SolverCallback & callback, + const ID & solver_id = ""); + + /// Initialize a time solver that can be used afterwards with its id void getNewSolver(const ID & solver_id, TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type = _nls_auto); /// set an integration scheme for a given dof and a given solver void setIntegrationScheme(const ID & solver_id, const ID & dof_id, const IntegrationSchemeType & integration_scheme_type, IntegrationScheme::SolutionType solution_type = IntegrationScheme::_not_defined); /// set an externally instantiated integration scheme void setIntegrationScheme(const ID & solver_id, const ID & dof_id, IntegrationScheme & integration_scheme, IntegrationScheme::SolutionType solution_type = IntegrationScheme::_not_defined); /* ------------------------------------------------------------------------ */ /* SolverCallback interface */ /* ------------------------------------------------------------------------ */ public: /// Predictor interface for the callback void predictor() override; /// Corrector interface for the callback void corrector() override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Default time step solver to instantiate for this model virtual TimeStepSolverType getDefaultSolverType() const; /// Default configurations for a given time step solver virtual ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /// get access to the internal dof manager DOFManager & getDOFManager() { return *this->dof_manager; } /// get the time step of a given solver Real getTimeStep(const ID & solver_id = "") const; /// set the time step of a given solver virtual void setTimeStep(Real time_step, const ID & solver_id = ""); /// set the parameter 'param' of the solver 'solver_id' // template // void set(const ID & param, const T & value, const ID & solver_id = ""); /// get the parameter 'param' of the solver 'solver_id' // const Parameter & get(const ID & param, const ID & solver_id = "") const; /// answer to the question "does the solver exists ?" bool hasSolver(const ID & solver_id) const; /// changes the current default solver void setDefaultSolver(const ID & solver_id); /// is a default solver defined bool hasDefaultSolver() const; /// is an integration scheme set for a given solver and a given dof bool hasIntegrationScheme(const ID & solver_id, const ID & dof_id) const; TimeStepSolver & getTimeStepSolver(const ID & solver_id = ""); NonLinearSolver & getNonLinearSolver(const ID & solver_id = ""); const TimeStepSolver & getTimeStepSolver(const ID & solver_id = "") const; const NonLinearSolver & getNonLinearSolver(const ID & solver_id = "") const; private: TimeStepSolver & getSolver(const ID & solver_id); const TimeStepSolver & getSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: ModelType model_type; private: ID parent_id; UInt parent_memory_id; /// Underlying mesh Mesh & mesh; /// Underlying dof_manager (the brain...) std::unique_ptr dof_manager; /// Default time step solver to use ID default_solver_id; }; struct ModelSolverOptions { NonLinearSolverType non_linear_solver_type; std::map integration_scheme_type; std::map solution_type; }; } // akantu #endif /* __AKANTU_MODEL_SOLVER_HH__ */