diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index a851126cf..ad6cc8aea 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,498 +1,429 @@ /** * @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) : Memory(id, memory_id), Parsable(ParserType::_contact_detector, id), mesh(mesh) { AKANTU_DEBUG_IN(); this->spatial_dimension = mesh.getSpatialDimension(); this->positions = positions; this->mesh.fillNodesToElements(this->spatial_dimension - 1); this->parseSection(); this->computeMaximalDetectionDistance(); 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"); surfaces[Surface::slave ] = section.getParameterValue("slave"); } /* -------------------------------------------------------------------------- */ void ContactDetector::search(std::map & contact_map) { SpatialGrid master_grid(spatial_dimension); SpatialGrid slave_grid(spatial_dimension); this->globalSearch(slave_grid, master_grid); this->localSearch(slave_grid, master_grid); this->constructContactMap(contact_map); } /* -------------------------------------------------------------------------- */ void ContactDetector::globalSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid) { 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); AKANTU_DEBUG_INFO( "Intersection BBox " << bbox_intersection ); Vector center(spatial_dimension); bbox_intersection.getCenter(center); Vector spacing(spatial_dimension); this->computeCellSpacing(spacing); master_grid.setCenter(center); master_grid.setSpacing(spacing); this->constructGrid(master_grid, bbox_intersection, master_list); slave_grid.setCenter(center); slave_grid.setSpacing(spacing); this->constructGrid(slave_grid, bbox_intersection, slave_list); // 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 elements 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. } /* -------------------------------------------------------------------------- */ void ContactDetector::localSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid) { // 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 contact_pairs.clear(); /// find the closet master node for each slave node for (auto && cell_id : slave_grid) { /// loop over all the slave nodes of the current cell for (auto && slave_node: slave_grid.getCell(cell_id)) { bool pair_exists = false; Vector pos(spatial_dimension); for (UInt s: arange(spatial_dimension)) pos(s) = this->positions(slave_node, s); Real closet_distance = std::numeric_limits::max(); UInt closet_master_node; /// loop over all the neighboring cells of the current cell for (auto && neighbor_cell : cell_id.neighbors()) { /// loop over the data of neighboring cells from master grid for (auto && master_node : master_grid.getCell(neighbor_cell)) { Vector pos2(spatial_dimension); for (UInt s: arange(spatial_dimension)) pos2(s) = this->positions(master_node, s); Real distance = pos.distance(pos2); if (distance <= closet_distance) { closet_master_node = master_node; closet_distance = distance; pair_exists = true; } } } if (pair_exists) { contact_pairs.push_back( std::make_pair(slave_node, closet_master_node)); } } } } /* -------------------------------------------------------------------------- */ void ContactDetector::constructContactMap(std::map & contact_map) { + auto surface_dimension = spatial_dimension - 1; + auto get_index = [&](auto & gaps, auto & projections) { UInt index; Real gap_min = std::numeric_limits::max(); UInt counter = 0; for (auto && values : zip(gaps, - make_view(projections, spatial_dimension -1))) { + make_view(projections, surface_dimension))) { auto & gap = std::get<0>(values); auto & projection = std::get<1>(values); - + + /// todo adhoc fix to assign a master element in case the + /// projection does not lie in the extended element. As it is + /// tolerance based bool is_valid = this->checkValidityOfProjection(projection); if (is_valid and gap <= gap_min) { gap_min = gap; index = counter; } counter++; } - /// todo adhoc fix to assign a master element in case the - /// projection does not lie in the extended element. As it is - /// tolerance based if (index >= gaps.size()) { auto gap_min_it = std::min_element(gaps.begin(), gaps.end()); auto index_it = std::find(gaps.begin(), gaps.end(), *gap_min_it); index = *index_it; } return index; }; auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn = this->mesh.getConnectivity(master); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i-1]; } return elem_conn; }; - std::cerr << contact_pairs.size() << std::endl; for (auto & pairs : contact_pairs) { const auto & slave_node = pairs.first; const auto & master_node = pairs.second; Array elements; this->mesh.getAssociatedElements(master_node, elements); - Array gaps(elements.size(), 1, "gaps"); - Array normals(elements.size(), spatial_dimension, "normals"); - Array projections(elements.size(), spatial_dimension -1, "projections"); + Array gaps(elements.size(), 1, "gaps"); + Array normals(elements.size(), spatial_dimension, "normals"); + Array projections(elements.size(), surface_dimension, "projections"); this->computeOrthogonalProjection(slave_node, elements, normals, gaps, projections); auto index = get_index(gaps, projections); auto connectivity = get_connectivity(slave_node, elements[index]); contact_map[slave_node].setMaster(elements[index]); contact_map[slave_node].setGap(gaps[index]); contact_map[slave_node].setNormal(Vector(normals.begin(spatial_dimension)[index], true)); - contact_map[slave_node].setProjection(Vector(projections.begin(spatial_dimension - 1)[index], true)); + contact_map[slave_node].setProjection(Vector(projections.begin(surface_dimension)[index], true)); contact_map[slave_node].setConnectivity(connectivity); - - /// number of surface tangents will be equal to dimension of - /// surface i.e. spatial_dimension - 1 and nb of components will - /// still be equal to spatial dimension - Matrix tangents(spatial_dimension - 1, spatial_dimension); + + Matrix tangents(surface_dimension, spatial_dimension); this->computeTangentsOnElement(contact_map[slave_node].master, contact_map[slave_node].projection, tangents); contact_map[slave_node].setTangent(tangents); } - + + contact_pairs.clear(); } /* -------------------------------------------------------------------------- */ void ContactDetector::computeOrthogonalProjection(const UInt & node, const Array & elements, Array & normals, Array & gaps, Array & projections) { Vector query(spatial_dimension); for (UInt s: arange(spatial_dimension)) query(s) = this->positions(node, s); for (auto && values : zip( elements, gaps, make_view(normals , spatial_dimension), make_view(projections, spatial_dimension - 1))) { const auto & element = std::get<0>(values); auto & gap = std::get<1>(values); auto & normal = std::get<2>(values); auto & projection = std::get<3>(values); this->computeNormalOnElement(element, normal); Vector real_projection(spatial_dimension); this->computeProjectionOnElement(element, normal, query, 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 = distance.dot(normal); Real tolerance = 1e-8; - + + // todo adhoc fix does not work always + normal *= -1.0; + /// todo: adhoc fix to ensure that normal is always into the slave /// surface. However, it doesnot work if gap is 0 as cos angle is /// a nan value //if (std::abs(cos_angle + 1) <= tolerance) { // normal *= -1.0; //} if (std::abs(cos_angle - 1) <= tolerance and detection_type == _explicit) { gap *= -1; } } } -/* -------------------------------------------------------------------------- */ -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)); Real alpha = (query - point).dot(normal); real_projection = query - alpha * normal; this->computeNaturalProjection(element, real_projection, 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, 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 } -/* -------------------------------------------------------------------------- */ -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::computeTangentsOnElement(const Element & el, Vector & projection, Matrix & tangents) { const ElementType & type = el.type; UInt nb_nodes_master = Mesh::getNbNodesPerElement(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 coords(spatial_dimension, nb_nodes_master); coordinatesOfElement(el, coords); tangents.mul(shapes_derivatives, coords); auto temp_tangents = tangents.transpose(); for (UInt i = 0; i < spatial_dimension -1; ++i) { auto temp = Vector(temp_tangents(i)); temp_tangents(i) = temp.normalize(); } tangents = temp_tangents.transpose(); } - -/* -------------------------------------------------------------------------- */ -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 e0bd15ce7..494db8d35 100644 --- a/src/model/contact_mechanics/contact_detector.hh +++ b/src/model/contact_mechanics/contact_detector.hh @@ -1,220 +1,203 @@ /** * @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 &, const ID & id = "contact_detector", UInt memory_id = 0); - ContactDetector(Mesh &, Array positions, - const ID & id = "contact_detector", + ContactDetector(Mesh &, Array positions, const ID & id = "contact_detector", UInt memory_id = 0); ~ContactDetector() = default; /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ public: - /// performs the all search steps + /// performs all search steps void search(std::map &); /// performs global spatial search to construct spatial grids void globalSearch(SpatialGrid &, SpatialGrid &); /// performs local search to find closet master node to a slave node void localSearch(SpatialGrid &, SpatialGrid &); /// constructs contact map for a given pair of slave and master node void constructContactMap(std::map &); private: /// reads the input file to get contact detection options void parseSection(); - /// extracts vectors which forms the plane of element - void vectorsAlongElement(const Element & /* element id */, - Matrix & /* vectors matrix */); - - /// computes orthogonal projection on master elements - void computeOrthogonalProjection(const UInt & /* slave node */, - const Array & /* master elements */, - Array & /* normals */, - Array & /* gaps */, - Array & /* projections */); + /// computes the orthogonal projection on master elements + void computeOrthogonalProjection(const UInt & , const Array &, + Array &, Array &, Array &); - /// computes normal on an element - void computeNormalOnElement(const Element & /* element id */, - Vector & /* normal vector */); - /// computes tangents on a given natural coordinate - void computeTangentsOnElement(const Element &, Vector &, - Matrix &); + void computeTangentsOnElement(const Element &, Vector &, Matrix &); /// 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 */); + void computeProjectionOnElement(const Element &, const Vector &, + const Vector &, Vector &, Vector &); /// computes natural projection of a real projection - void computeNaturalProjection(const Element & /* element */, - Vector & /* real projection */, - Vector & /* natural projection */); - - /// compute normal projection of slave coord on a given element - void normalProjection(const Element & el, const Vector & slave_coord, - Vector & natural_coord, Real & tolerance); + void computeNaturalProjection(const Element &, Vector &, Vector &); /* ------------------------------------------------------------------------ */ /* Inline Methods */ /* ------------------------------------------------------------------------ */ public: /// checks whether the natural projection is valid or not inline bool checkValidityOfProjection(Vector & ); /// extracts the coordinates of an element inline void coordinatesOfElement(const Element & , Matrix & ); /// computes the optimal cell size for grid inline void computeCellSpacing(Vector & ); /// constructs a grid containing nodes lying within bounding box inline void constructGrid(SpatialGrid &, BBox &, const Array &); /// constructs the bounding box based on nodes list inline void constructBoundingBox(BBox &, const Array &); /// get the surface id template inline std::string getSurfaceId(); /// set the surface id template inline void setSurfaceId(const std::string); /// computes the maximum in radius for a given mesh inline void computeMaximalDetectionDistance(); /// constructs the connectivity for a contact element inline Vector constructConnectivity(UInt &, const Element &); + /// computes normal on an element + inline void computeNormalOnElement(const Element &, Vector & ); + + /// extracts vectors which forms the plane of element + inline void vectorsAlongElement(const Element &, Matrix & ); + + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// returns the maximum detection distance AKANTU_GET_MACRO(MaximumDetectionDistance, max_dd, Real); /// sets the maximum detection distance AKANTU_SET_MACRO(MaximumDetectionDistance, max_dd, Real); /// returns the bounding box extension AKANTU_GET_MACRO(MaximumBoundingBox, max_bb, Real); /// sets the bounding box extension AKANTU_SET_MACRO(MaximumBoundingBox, max_bb, Real); /// returns the positions AKANTU_GET_MACRO(Positions, positions, Array); /// sets the positions AKANTU_SET_MACRO(Positions, positions, Array); /* ------------------------------------------------------------------------ */ /* 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; /// contact pair slave node to closet master node std::vector > contact_pairs; /// contains the updated positions of the nodes Array positions; /// type of detection explicit/implicit DetectionType detection_type; }; } // namespace akantu #include "contact_detector_inline_impl.cc" #endif /* __AKANTU_CONTACT_DETECTOR_HH__ */ diff --git a/src/model/contact_mechanics/contact_detector_inline_impl.cc b/src/model/contact_mechanics/contact_detector_inline_impl.cc index 3b4ed2493..ae4f019ce 100644 --- a/src/model/contact_mechanics/contact_detector_inline_impl.cc +++ b/src/model/contact_mechanics/contact_detector_inline_impl.cc @@ -1,191 +1,233 @@ /** * @file contact_detection.hh * * @author Mohit Pundir * * @date creation: Mon Apr 29 2019 * @date last modification: Mon Apr 29 2019 * * @brief inine implementation of the contact detector 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 "contact_detector.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_DETECTOR_INLINE_IMPL_CC__ #define __AKANTU_CONTACT_DETECTOR_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline bool ContactDetector::checkValidityOfProjection(Vector & projection) { UInt nb_xi_inside = 0; Real epsilon = 1e-3; for (auto xi : projection) { if (xi >= -1.0 - epsilon and xi <= 1.0 + epsilon) nb_xi_inside++; } if (nb_xi_inside == projection.size()) return true; return false; } /* -------------------------------------------------------------------------- */ inline 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); } } } /* -------------------------------------------------------------------------- */ inline void ContactDetector::computeCellSpacing(Vector & spacing) { for (UInt s: arange(spatial_dimension)) spacing(s) = std::sqrt(2.0) * max_dd; } /* -------------------------------------------------------------------------- */ inline 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); 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; } - - std::cerr << bbox << std::endl; AKANTU_DEBUG_INFO("BBox" << bbox); } /* -------------------------------------------------------------------------- */ inline 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); } /* -------------------------------------------------------------------------- */ template inline std::string ContactDetector::getSurfaceId() { return surfaces[id]; } /* -------------------------------------------------------------------------- */ template inline void ContactDetector::setSurfaceId(const std::string name) { surfaces[id] = name; } /* -------------------------------------------------------------------------- */ inline void ContactDetector::computeMaximalDetectionDistance() { AKANTU_DEBUG_IN(); Real elem_size; Real max_elem_size = std::numeric_limits::min(); auto & master_group = mesh.getElementGroup(surfaces[Surface::master]); 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); 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); elem_size = FEEngine::getElementInradius(elem_coords, type); max_elem_size = std::max(max_elem_size, elem_size); } AKANTU_DEBUG_INFO("The maximum element size : " << max_elem_size ); } this->max_dd = max_elem_size; this->max_bb = max_elem_size; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline Vector ContactDetector::constructConnectivity(UInt & slave, const Element & master) { Vector master_conn = this->mesh.getConnectivity(master); Vector elem_conn(master_conn.size() + 1); elem_conn[0] = slave; for (UInt i = 1; i < elem_conn.size(); ++i) { elem_conn[i] = master_conn[i-1]; } return elem_conn; } + +/* -------------------------------------------------------------------------- */ +inline 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); } + } + +} + +/* -------------------------------------------------------------------------- */ +inline 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); } + } + +} + /* -------------------------------------------------------------------------- */ } //akantu #endif /* __AKANTU_CONTACT_DETECTOR_INLINE_IMPL_CC__ */ diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 67b598811..557a27e51 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,645 +1,599 @@ /** * @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 "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, Array & positions, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, 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->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, current_positions, id + ":contact_detector"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : ContactMechanicsModel(mesh, mesh.getNodes(), dim, id, memory_id, dof_manager, model_type) {} /* -------------------------------------------------------------------------- */ 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) { auto & dof_manager = this->getDOFManager(); /* -------------------------------------------------------------------------- */ // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->normals, spatial_dimension, "normals"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->nodal_area, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "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(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the contact forces"); this->internal_force->clear(); // assemble the forces due to contact 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); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { this->contact_map.clear(); - this->detector->search(this->contact_map); - this->assembleFieldsFromContactMap(); - this->computeNodalAreas(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search(Array & increment) { this->contact_map.clear(); - this->detector->search(this->contact_map); for (auto & entry : contact_map) { auto & element = entry.second; const auto & connectivity = element.connectivity; const auto & normal = element.normal; auto slave_node = connectivity[0]; Vector u_slave(spatial_dimension); auto master_node = connectivity[1]; Vector u_master(spatial_dimension); for (UInt s : arange(spatial_dimension)) { u_slave(s) = increment(slave_node, s); u_master(s) = increment(master_node, s); } - /*auto u = (u_master.norm() > u_slave.norm()) ? u_master * -1.0 : u_slave; - - const auto & normal = element.normal; - Real uv = Math::vectorDot(u.storage(), normal.storage(), spatial_dimension); - - std::cerr << u << " " << normal << std::endl; - std::cerr << uv << " " << element.gap << " " << uv + element.gap << std::endl; - - if (uv + element.gap <= 0) { - element.gap = abs(uv + element.gap); - } else { - element.gap = 0.0; - }*/ - + // todo check this initial guess auto u = (u_master.norm() > u_slave.norm()) ? u_master : u_slave; - Real uv = Math::vectorDot(u.storage(), normal.storage(), spatial_dimension); - uv = abs(uv); - - std::cerr << uv << " " << element.gap << " " << element.gap - uv << std::endl; - + if (element.gap - uv <= 0) { element.gap = abs(element.gap - uv); - } else { - element.gap = 0.0; - } - - /*auto u = u_slave - u_master; - auto u_n = u.dot(normal); - u_n = abs(u_n); - - auto positions = detector->getPositions(); - Vector slave_pos(spatial_dimension); - Vector master_pos(spatial_dimension); - - for (auto s : arange(spatial_dimension)) { - slave_pos(s) = positions(slave_node, s); - master_pos(s) = positions(master_node, s); - } - - auto sign = (slave_pos - master_pos).dot(normal); - sign/= abs(sign); - - auto g_trial = element.gap - sign * u_n; - if (g_trial <= 0 or sign < 0) { - element.gap = abs(g_trial); } else { - element.gap = 0; - }*/ - - + element.gap = 0.0; + } } this->assembleFieldsFromContactMap(); this->computeNodalAreas(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleFieldsFromContactMap() { if (this->contact_map.empty()) AKANTU_ERROR( "Contact map is empty, Please run search before assembling the fields"); for (auto & entry : contact_map) { const auto & element = entry.second; auto connectivity = element.connectivity; auto node = connectivity(0); (*gaps)[node] = element.gap; for (UInt i = 0; i < spatial_dimension; ++i) (*normals)(node, i) = element.normal[i]; } } /* -------------------------------------------------------------------------- */ template void ContactMechanicsModel::computeNodalAreas() { this->nodal_area->clear(); this->external_force->clear(); this->applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), this->detector->getSurfaceId()); for (auto && tuple : zip(*nodal_area, make_view(*external_force, spatial_dimension))) { auto & area = std::get<0>(tuple); auto & force = std::get<1>(tuple); for (auto & f : force) area += pow(f, 2); area = sqrt(area); } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ 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) { 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")); } 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*/, bool /*padding_flag*/) { dumper::Field * field = nullptr; 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->internal_force; real_nodal_fields["blocked_dofs"] = this->blocked_dofs; real_nodal_fields["normals"] = this->normals; real_nodal_fields["gaps"] = this->gaps; real_nodal_fields["areas"] = this->nodal_area; 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; 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/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc index dd8979571..3a97e9f90 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc @@ -1,321 +1,348 @@ /** * @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, "Normal penalty parameter"); this->registerParam("epsilon_t", epsilon_t, Real(0.), _pat_parsable | _pat_modifiable, "Tangential penalty parameter"); } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeNormalForce(Vector & force, Vector & n, ContactElement & element) { force.clear(); Real tn = element.gap * epsilon; tn = macaulay(tn); for (UInt i : arange(force.size())) { force[i] += n[i] * tn; } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeFrictionalForce(Vector & force, Array & d_alpha, ContactElement & element) { Matrix m_alpha_beta(spatial_dimension - 1, spatial_dimension - 1); ResolutionUtils::computeMetricTensor(m_alpha_beta, element.tangents); computeFrictionalTraction(m_alpha_beta, element); auto & traction = element.traction; for (auto && values: zip(traction, 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::computeNormalModuli(Matrix & ke, Array & n_alpha, Array & d_alpha, Vector & n, ContactElement & element) { Real tn = element.gap * epsilon; tn = macaulay(tn); Matrix n_mat(n.storage(), n.size(), 1); - ke.mul(n_mat, n_mat); ke *= epsilon * heaviside(element.gap); 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); + tmp2.mul(ds_mat, ns_mat); ke -= (tmp1 + tmp2) * tn; } + + auto surface_dimension = spatial_dimension - 1; + Matrix m_alpha_beta(surface_dimension, surface_dimension); + ResolutionUtils::computeMetricTensor(m_alpha_beta, element.tangents); + + for (auto && values : + zip(arange(surface_dimension), + make_view(n_alpha, n_alpha.size()))) { + + auto & s = std::get<0>(values); + auto & n_s = std::get<1>(values); + Matrix ns_mat(n_s.storage(), n_s.size(), 1); + + for (auto && tuple : + zip(arange(surface_dimension), + make_view(n_alpha, n_alpha.size()))) { + auto & t = std::get<0>(tuple); + auto & n_t = std::get<1>(tuple); + + Matrix nt_mat(n_t.storage(), n_t.size(), 1); + + Matrix tmp3(n_s.size(), n_s.size()); + tmp3.mul(ns_mat, nt_mat); + + tmp3 *= m_alpha_beta(s, t) * tn * element.gap; + ke += tmp3; + } + } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeFrictionalModuli(Matrix & /*ke*/, Array & t_alpha_beta, Array & n_alpha_beta, Array & /*n_alpha*/, Array & d_alpha, Matrix & phi, Vector & n, ContactElement & element) { auto k_common = computeCommonModuli(t_alpha_beta, n_alpha_beta, d_alpha, n, element); const auto & type = element.master.type; const auto & conn = element.connectivity; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; Matrix m_alpha_beta(surface_dimension, surface_dimension); ResolutionUtils::computeMetricTensor(m_alpha_beta, element.tangents); Array g_alpha(conn.size() * spatial_dimension, surface_dimension); ResolutionUtils::computeGalpha(g_alpha, t_alpha_beta, d_alpha, phi, element); /*Matrix k_t; bool stick = computeFrictionalTraction(m_alpha_beta, element); if(stick) k_t = computeStickModuli(g_alpha, d_alpha, m_alpha_beta); else k_t = computeSlipModuli(g_alpha, d_alpha, m_alpha_beta, element);*/ } /* -------------------------------------------------------------------------- */ bool ResolutionPenalty::computeFrictionalTraction(Matrix& m_alpha_beta, ContactElement & element) { Real tn = element.gap * epsilon; tn = macaulay(tn); auto delta_xi = element.projection - element.previous_projection; Vector trial_traction; trial_traction.mul(m_alpha_beta, delta_xi, epsilon); trial_traction += element.traction; auto trial_slip_function = trial_traction.norm() - mu * tn; bool stick = false; if (trial_slip_function <= 0) { element.traction = trial_traction; stick = true; } else{ element.traction = mu * tn * trial_traction / trial_traction.norm(); } return stick; } /* -------------------------------------------------------------------------- */ Array ResolutionPenalty::computeCommonModuli(Array & t_alpha_beta, Array & n_alpha_beta, Array & d_alpha, Vector & n, ContactElement & element) { Array kt_alpha(spatial_dimension -1, d_alpha.size() * d_alpha.size(), "k_T_alpha"); //auto t_alpha_beta_size = t_alpha_beta.size() * (spatial_dimension - 1); //auto & tangents = element.tangents; /*for(auto && values : zip(tangents.transpose(), make_view(kt_alpha, kt_alpha.size()), make_view(t_alpha_beta, t_alpha_beta_size), make_view(n_alpha_beta, n_alpha_beta_size))) { auto & tangent_s = std::get<0>(values); auto & kt_s = std::get<1>(values); auto & t_alpha_s = std::get<2>(values); auto & n_alpha_s = std::get<3>(values); Matrix kt_s_mat(kt_s.storage(), d_alpha.size(), d_alpha.size()); // loop over beta for(auto && tuple : zip(make_view(d_alpha, d_alpha.size()), make_view(n_alpha_ ))) { auto & d_s = std::get<0>(tuple); Matrix tmp(d_s.size(), d_s.size()); // loop over gamma for(auto && entry : make_view(d_alpha, d_alpha.size())) { auto & d_t = std::get<0>(entry); // compute constant Matrix tmp2(d_t.size(), d_t.size()); tmp2.mul(d_s, d_t); kt_s_mat += tmp2; } } }*/ return kt_alpha; } /* -------------------------------------------------------------------------- */ Matrix ResolutionPenalty::computeStickModuli(Array & g_alpha, Array & d_alpha, Matrix & m_alpha_beta) { Matrix k_stick(d_alpha.size(), d_alpha.size()); for (auto && values : zip(arange(d_alpha.getNbComponent()), make_view(d_alpha, d_alpha.size()), make_view(g_alpha, g_alpha.size()))) { auto & s = std::get<0>(values); auto & d_s = std::get<1>(values); auto & g_s = std::get<2>(values); Matrix ds_mat(d_s.storage(), d_s.size(), 1); Matrix gs_mat(g_s.storage(), g_s.size(), 1); Matrix tmp1(d_s.size(), d_s.size()); tmp1.mul(ds_mat, gs_mat); k_stick += tmp1; for(auto && tuple : enumerate(make_view(d_alpha, d_alpha.size()))) { auto & t = std::get<0>(tuple); auto & d_t = std::get<1>(tuple); Matrix dt_mat(d_t.storage(), d_t.size(), 1); Matrix tmp2(d_t.size(), d_t.size()); tmp2.mul(ds_mat, dt_mat); k_stick += tmp2 * m_alpha_beta(s, t); } } k_stick *= epsilon_t; return k_stick; } /* -------------------------------------------------------------------------- */ Matrix ResolutionPenalty::computeSlipModuli(Array & g_alpha, Array & d_alpha, Matrix & m_alpha_beta, ContactElement & element) { Real tn = element.gap * epsilon; tn = macaulay(tn); Real factor; factor = epsilon_t * mu * tn; auto p_t = element.traction; p_t /= p_t.norm(); Matrix k_slip(d_alpha.size(), d_alpha.size()); /* // loop for alpha for(auto && value : make_view(d_alpha, d_alpha.size())) { auto & d_s = std::get<0>(value); // loop for beta for(auto && tuple : zip(arange(spatial_dimension - 1), make_view(d_alpha, d_alpha.size()), make_view(g_alpha, g_alpha.size()))) { auto & beta = std::get<0>(tuple); auto & d_beta = std::get<1>(tuple); auto & g_beta = std::get<2>(tuple); // loop for gamma for(auto && entry : zip(arange(spatial_dimension - 1), make_view(d_alpha, d_alpha.size()))) { auto & gamma = std::get<0>(entry); auto & d_gamma = std::get<1>(entry); } } }*/ return k_slip; } INSTANTIATE_RESOLUTION(penalty, ResolutionPenalty); } // akantu diff --git a/src/model/model_couplers/coupler_solid_contact.cc b/src/model/model_couplers/coupler_solid_contact.cc index a071bff4c..f5e47fd77 100644 --- a/src/model/model_couplers/coupler_solid_contact.cc +++ b/src/model/model_couplers/coupler_solid_contact.cc @@ -1,411 +1,412 @@ /** * @file coupler_solid_contact_explicit.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 * @date last modification: Thu May 22 2019 * * @brief class for coupling of solid mechanics and conatct 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 "coupler_solid_contact.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidContact::CouplerSolidContact(Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("CouplerSolidContact", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModel(mesh, Model::spatial_dimension, "solid_mechanics_model", 0, this->dof_manager); contact = new ContactMechanicsModel(mesh, Model::spatial_dimension, "contact_mechanics_model", 0, this->dof_manager); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidContact::~CouplerSolidContact() {} /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidContact::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::initSolver(TimeStepSolverType, NonLinearSolverType) { DOFManager & dof_manager = this->getDOFManager(); } /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidContact::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 CouplerSolidContact::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 CouplerSolidContact::assembleResidual() { // computes the internal forces this->assembleInternalForces(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); auto & contact_map = contact->getContactMap(); for (auto & pair: contact_map) { auto & connectivity = pair.second.connectivity; for (auto node : connectivity) { for (auto s : arange(spatial_dimension)) { external_force(node, s) = contact_force(node, s); } } } /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", - internal_force, 1); + internal_force, 1); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); auto & contact_map = contact->getContactMap(); for (auto & pair: contact_map) { auto & connectivity = pair.second.connectivity; for (auto node : connectivity) { for (auto s : arange(spatial_dimension)) { external_force(node, s) = contact_force(node, s); } } } if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("displacement", internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::beforeSolveStep() { Array increment(0, Model::spatial_dimension); auto & previous_displacement = solid->getPreviousDisplacement(); auto & displacement = solid->getDisplacement(); auto disp_it = displacement.begin(Model::spatial_dimension); auto disp_end = displacement.end(Model::spatial_dimension); auto prev_it = previous_displacement.begin(Model::spatial_dimension); for (; disp_it != disp_end; ++disp_it, ++prev_it) { increment.push_back(*disp_it - *prev_it); } contact->search(displacement); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::afterSolveStep() { Array current_positions(0, Model::spatial_dimension); auto & positions = mesh.getNodes(); auto & displacement = solid->getDisplacement(); auto pos_it = positions.begin(Model::spatial_dimension); auto pos_end = positions.end(Model::spatial_dimension); auto disp_it = displacement.begin(Model::spatial_dimension); for (; pos_it != pos_end; ++pos_it, ++disp_it) { current_positions.push_back(*pos_it + *disp_it); } contact->setPositions(current_positions); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::predictor() {} /* -------------------------------------------------------------------------- */ void CouplerSolidContact::corrector() {} /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidContact::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleLumpedMatrix(const ID & /*matrix_id*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); contact->assembleInternalForces(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); contact->assembleStiffnessMatrix(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ dumper::Field * CouplerSolidContact::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { dumper::Field * field = nullptr; field = solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * CouplerSolidContact::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { dumper::Field * field = nullptr; if (field_name == "contact_force" or field_name == "normals" or field_name == "gaps") field = contact->createNodalFieldReal(field_name, group_name, padding_flag); else field = solid->createNodalFieldReal(field_name, group_name, padding_flag); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * CouplerSolidContact::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { dumper::Field * field = nullptr; field = solid->createNodalFieldBool(field_name, group_name, padding_flag); return field; } #else /* -------------------------------------------------------------------------- */ dumper::Field * CouplerSolidContact::createElementalField(const std::string &, const std::string &, bool, const UInt &, const ElementKind &) { return nullptr; } /* ----------------------------------------------------------------------- */ dumper::Field * CouplerSolidContact::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ dumper::Field * CouplerSolidContact::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void CouplerSolidContact::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidContact::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu diff --git a/src/model/non_linear_solver_newton_raphson.cc b/src/model/non_linear_solver_newton_raphson.cc index 67618fff6..c1fd7d232 100644 --- a/src/model/non_linear_solver_newton_raphson.cc +++ b/src/model/non_linear_solver_newton_raphson.cc @@ -1,192 +1,193 @@ /** * @file non_linear_solver_newton_raphson.cc * * @author Nicolas Richart * * @date creation: Tue Sep 15 2015 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of the default NonLinearSolver * * @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 "non_linear_solver_newton_raphson.hh" #include "communicator.hh" #include "dof_manager_default.hh" #include "solver_callback.hh" #include "sparse_solver_mumps.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NonLinearSolverNewtonRaphson::NonLinearSolverNewtonRaphson( DOFManagerDefault & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id, UInt memory_id) : NonLinearSolver(dof_manager, non_linear_solver_type, id, memory_id), dof_manager(dof_manager), solver(std::make_unique( dof_manager, "J", id + ":sparse_solver", memory_id)) { this->supported_type.insert(_nls_newton_raphson_modified); this->supported_type.insert(_nls_newton_raphson); this->supported_type.insert(_nls_linear); this->checkIfTypeIsSupported(); this->registerParam("threshold", convergence_criteria, 1e-10, _pat_parsmod, "Threshold to consider results as converged"); this->registerParam("convergence_type", convergence_criteria_type, _scc_solution, _pat_parsmod, "Type of convergence criteria"); this->registerParam("max_iterations", max_iterations, 10, _pat_parsmod, "Max number of iterations"); this->registerParam("error", error, _pat_readable, "Last reached error"); this->registerParam("nb_iterations", n_iter, _pat_readable, "Last reached number of iterations"); this->registerParam("converged", converged, _pat_readable, "Did last solve converged"); this->registerParam("force_linear_recompute", force_linear_recompute, true, _pat_modifiable, "Force reassembly of the jacobian matrix"); } /* -------------------------------------------------------------------------- */ NonLinearSolverNewtonRaphson::~NonLinearSolverNewtonRaphson() = default; /* ------------------------------------------------------------------------ */ void NonLinearSolverNewtonRaphson::solve(SolverCallback & solver_callback) { this->dof_manager.updateGlobalBlockedDofs(); solver_callback.predictor(); if (non_linear_solver_type == _nls_linear and solver_callback.canSplitResidual()) solver_callback.assembleMatrix("K"); this->assembleResidual(solver_callback); if (this->non_linear_solver_type == _nls_newton_raphson_modified || (this->non_linear_solver_type == _nls_linear && this->force_linear_recompute)) { solver_callback.assembleMatrix("J"); this->force_linear_recompute = false; } this->n_iter = 0; this->converged = false; + //this->dump(); if (this->convergence_criteria_type == _scc_residual) { this->converged = this->testConvergence(this->dof_manager.getResidual()); if (this->converged) return; } do { if (this->non_linear_solver_type == _nls_newton_raphson) solver_callback.assembleMatrix("J"); this->solver->solve(); solver_callback.corrector(); // EventManager::sendEvent(NonLinearSolver::AfterSparseSolve(method)); if (this->convergence_criteria_type == _scc_residual) { this->assembleResidual(solver_callback); this->converged = this->testConvergence(this->dof_manager.getResidual()); } else { this->converged = this->testConvergence(this->dof_manager.getGlobalSolution()); } if (this->convergence_criteria_type == _scc_solution and not this->converged) this->assembleResidual(solver_callback); - // this->dump(); + //this->dump(); this->n_iter++; AKANTU_DEBUG_INFO( "[" << this->convergence_criteria_type << "] Convergence iteration " << std::setw(std::log10(this->max_iterations)) << this->n_iter << ": error " << this->error << (this->converged ? " < " : " > ") << this->convergence_criteria); } while (not this->converged and this->n_iter < this->max_iterations); // this makes sure that you have correct strains and stresses after the // solveStep function (e.g., for dumping) if (this->convergence_criteria_type == _scc_solution) this->assembleResidual(solver_callback); if (this->converged) { // this->sendEvent(NonLinearSolver::ConvergedEvent(method)); } else if (this->n_iter == this->max_iterations) { AKANTU_CUSTOM_EXCEPTION(debug::NLSNotConvergedException( this->convergence_criteria, this->n_iter, this->error)); AKANTU_DEBUG_WARNING("[" << this->convergence_criteria_type << "] Convergence not reached after " << std::setw(std::log10(this->max_iterations)) << this->n_iter << " iteration" << (this->n_iter == 1 ? "" : "s") << "!"); } return; } /* -------------------------------------------------------------------------- */ bool NonLinearSolverNewtonRaphson::testConvergence(const Array & array) { AKANTU_DEBUG_IN(); const Array & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); UInt nb_degree_of_freedoms = array.size(); auto arr_it = array.begin(); auto bld_it = blocked_dofs.begin(); Real norm = 0.; for (UInt n = 0; n < nb_degree_of_freedoms; ++n, ++arr_it, ++bld_it) { bool is_local_node = this->dof_manager.isLocalOrMasterDOF(n); if ((!*bld_it) && is_local_node) { norm += *arr_it * *arr_it; } } dof_manager.getCommunicator().allReduce(norm, SynchronizerOperation::_sum); norm = std::sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something went wrong in the solve phase"); this->error = norm; return (error < this->convergence_criteria); } /* -------------------------------------------------------------------------- */ } // akantu diff --git a/test/test_model/test_contact_mechanics_model/CMakeLists.txt b/test/test_model/test_contact_mechanics_model/CMakeLists.txt index 570b5903d..1db816582 100644 --- a/test/test_model/test_contact_mechanics_model/CMakeLists.txt +++ b/test/test_model/test_contact_mechanics_model/CMakeLists.txt @@ -1,55 +1,60 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Tue Dec 18 2018 # @date last modification: Tue Dec 18 2018 # # @brief configuration for ContactMechanicsModel tests # # @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 . # # @section DESCRIPTION # #=============================================================================== add_akantu_test(test_detector "detection_test") add_akantu_test(test_resolution "resolution_test") add_akantu_test(test_coupler "coupler_test") -add_mesh (contact_2d contact_2d.geo 2 1) -add_mesh (contact_3d contact_3d.geo 3 1) +add_mesh (hertz_2d hertz_2d.geo 2 1) +add_mesh (hertz_3d hertz_3d.geo 3 1) add_mesh (circles_2d circles_2d.geo 2 1) +add_mesh (flat_2d flat_2d.geo 2 1) - -register_test(test_contact_mechanics - SOURCES test_contact_mechanics_model.cc - DEPENDS contact_3d +register_test(test_hertz_2d + SOURCES test_hertz_2d.cc + DEPENDS hertz_2d FILES_TO_COPY material.dat PACKAGE contact_mechanics ) -register_test(test_contact_mechanics_2d - SOURCES test_contact_mechanics_model_2d.cc - DEPENDS contact_2d +register_test(test_circles_2d + SOURCES test_circles_2d.cc + DEPENDS circles_2d FILES_TO_COPY material.dat PACKAGE contact_mechanics ) - -register_test(test_contact_mechanics_circles - SOURCES test_contact_mechanics_model.cc - DEPENDS circles_2d +register_test(test_hertz_3d + SOURCES test_hertz_3d.cc + DEPENDS hertz_3d FILES_TO_COPY material.dat PACKAGE contact_mechanics ) +register_test(test_flat_2d + SOURCES test_flat_2d.cc + DEPENDS flat_2d + FILES_TO_COPY material.dat + PACKAGE contact_mechanics + ) diff --git a/test/test_model/test_contact_mechanics_model/circles_2d.geo b/test/test_model/test_contact_mechanics_model/circles_2d.geo index f29094576..e648dc710 100644 --- a/test/test_model/test_contact_mechanics_model/circles_2d.geo +++ b/test/test_model/test_contact_mechanics_model/circles_2d.geo @@ -1,43 +1,43 @@ -cl1 = 0.0025; -cl2 = 0.0025; +cl1 = 0.001; +cl2 = 0.005; cl3 = 0.005; Dy = 0.0; radius = 0.1; y = 0.1; -epsilon = 0.0; -Point(1) = {0, y, 0, cl1}; -Point(2) = {radius, radius + y, 0, cl2}; -Point(3) = {-radius, radius + y, 0, cl2}; +epsilon = -1e-6; +Point(1) = {0, y + epsilon, 0, cl1}; +Point(2) = {radius, radius + y + epsilon, 0, cl2}; +Point(3) = {-radius, radius + y + epsilon, 0, cl2}; -Point(11) = {0, y, 0, cl1}; +Point(11) = {0, y, 0, cl2}; Point(12) = {radius, -radius + y, 0, cl2}; Point(13) = {-radius,- radius + y, 0, cl2}; Point(8) = {0, radius + y, 0, cl2}; Point(18) = {0, -radius + y, 0, cl2}; Circle(1) = {3, 8, 1}; Circle(2) = {1, 8, 2}; Circle(11) = {13, 18, 11}; Circle(12) = {11, 18, 12}; Line(3) = {2, 8}; Line(13) = {8, 3}; Line(30) = {12, 18}; Line(31) = {18, 13}; Line Loop(9) = {2, 3, 13, 1}; Plane Surface(9) = {9}; Line Loop(19) = {-12, -30, -31, -11}; Plane Surface(19) = {19}; -Physical Line("rigid") = {11, 12}; -Physical Line("contact_surface") = {1, 2}; +Physical Line("contact_bottom") = {11, 12}; +Physical Line("contact_top") = {1, 2}; Physical Line("bottom") = {30, 31}; Physical Line("top") = {3, 13}; Physical Surface("bot_body") = {19}; Physical Surface("top_body") = {9}; \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/flat_2d.geo b/test/test_model/test_contact_mechanics_model/flat_2d.geo new file mode 100644 index 000000000..5085f88dd --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/flat_2d.geo @@ -0,0 +1,33 @@ +a = 1.0; +r = a; +gap = 0; +nbpoints = 1.0; +cl1 = a/(nbpoints); +cl2 = r/(nbpoints); + +Point(1) = {-a/2, 0 + gap, 0, cl1}; +Point(2) = {a/2, 0 + gap, 0, cl1}; +Point(3) = {a/2, a/2., 0, cl1}; +Point(4) = {-a/2, a/2., 0, cl1}; +Point(5) = {r/2, 0, 0, cl2}; +Point(6) = {-r/2, 0, 0, cl2}; +Point(7) = {-r/2, -r/2, 0, 10*cl2}; +Point(8) = { r/2, -r/2, 0, 10*cl2}; +Line(2) = {1, 2}; +Line(3) = {2, 3}; +Line(4) = {3, 4}; +Line(5) = {4, 1}; +Line(6) = {7, 8}; +Line(7) = {8, 5}; +Line(8) = {5, 6}; +Line(9) = {6, 7}; +Line Loop(1) = {2, 3, 4, 5}; +Line Loop(2) = {6, 7, 8, 9}; +Plane Surface(1) = {1}; +Plane Surface(2) = {2}; +Physical Line("contact_top") = {2}; +Physical Line("contact_bottom") = {8}; +Physical Line("top") = {4}; +Physical Line("bottom") = {6}; +Physical Surface("bot_body") = {2}; +Physical Surface("top_body") = {1}; \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/contact_2d.geo b/test/test_model/test_contact_mechanics_model/hertz_2d.geo similarity index 85% rename from test/test_model/test_contact_mechanics_model/contact_2d.geo rename to test/test_model/test_contact_mechanics_model/hertz_2d.geo index d9bc359e6..81a055eef 100644 --- a/test/test_model/test_contact_mechanics_model/contact_2d.geo +++ b/test/test_model/test_contact_mechanics_model/hertz_2d.geo @@ -1,33 +1,33 @@ cl1 = 0.0025; cl2 = 0.0025; -cl3 = 0.005; +cl3 = 0.0025; Dy = 0.0; radius = 0.1; y = 0.1; epsilon = 0.0; Point(1) = {0, y, 0, cl1}; Point(2) = {radius, radius + y, 0, cl2}; Point(3) = {-radius, radius + y, 0, cl2}; Point(4) = {0.25, 0.1 -epsilon, 0, cl3}; Point(5) = {-0.25, 0.1-epsilon, 0, cl3}; Point(6) = {-0.25, -0.2, 0, 10*cl2}; Point(7) = {0.25, -0.2, 0, 10*cl2}; Point(8) = {0, radius + y, 0, cl2}; Circle(1) = {3, 8, 1}; Circle(2) = {1, 8, 2}; Line(3) = {2, 8}; Line(13) = {8, 3}; Line(4) = {6, 7}; Line(5) = {7, 4}; Line(6) = {4, 5}; Line(7) = {5, 6}; Line Loop(9) = {2, 3, 13, 1}; Plane Surface(9) = {9}; Line Loop(11) = {6, 7, 4, 5}; Plane Surface(11) = {11}; -Physical Line("flat") = {6}; -Physical Line("curved") = {1, 2}; +Physical Line("contact_bottom") = {6}; +Physical Line("contact_top") = {1, 2}; Physical Line("top") = {3, 13}; Physical Line("bottom") = {4}; Physical Surface("top_body") = {9}; Physical Surface("bot_body") = {11}; \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/contact_3d.geo b/test/test_model/test_contact_mechanics_model/hertz_3d.geo similarity index 83% rename from test/test_model/test_contact_mechanics_model/contact_3d.geo rename to test/test_model/test_contact_mechanics_model/hertz_3d.geo index 068722fb0..e03d1cc2f 100644 --- a/test/test_model/test_contact_mechanics_model/contact_3d.geo +++ b/test/test_model/test_contact_mechanics_model/hertz_3d.geo @@ -1,85 +1,85 @@ -cl1 = 0.005; -cl2 = 0.005; +cl1 = 0.01; +cl2 = 0.02; cl3 = 0.05; -cl4 = 0.005; -radius = 0.1; -depth = -0.1; +cl4 = 0.02; +radius = 0.25; +depth = -0.25; y = 0.0; Dz = 1; Point(1) = {0, 0, 0, cl1}; Point(2) = {radius/2, radius/2, 0, cl2}; Point(3) = {-radius/2, radius/2, 0, cl2}; Point(4) = {0, radius/2, 0, cl2}; Point(5) = {0, radius/2, -radius/2, cl2}; Point(6) = {0, radius/2, radius/2, cl2}; Point(12) = {0, y, 0, cl4}; Point(13) = {radius, y, radius, cl4}; Point(14) = {-radius, y, radius, cl4}; Point(15) = {-radius, y, -radius, cl4}; Point(16) = {radius, y, -radius, cl4}; -Point(17) = {radius, y + depth, radius, cl2}; -Point(18) = {-radius, y + depth, radius, cl2}; -Point(19) = {-radius, y + depth, -radius, cl2}; -Point(20) = {radius, y + depth, -radius, cl2}; +Point(17) = {radius, y + depth, radius, cl3}; +Point(18) = {-radius, y + depth, radius, cl3}; +Point(19) = {-radius, y + depth, -radius, cl3}; +Point(20) = {radius, y + depth, -radius, cl3}; Circle(1) = {3, 4, 1}; Circle(2) = {1, 4, 2}; Circle(3) = {6, 4, 1}; Circle(4) = {1, 4, 5}; Circle(5) = {3, 4, 6}; Circle(6) = {6, 4, 2}; Circle(7) = {2, 4, 5}; Circle(8) = {5, 4, 3}; Line Loop(1) = {3, 2, -6}; Ruled Surface(1) = {1}; Line Loop(2) = {-2, -7, 4}; Ruled Surface(2) = {2}; Line Loop(3) = {-4, -8, -1}; Ruled Surface(3) = {3}; Line Loop(4) = {1, -3, -5}; Ruled Surface(4) = {4}; Line Loop(5) = {6, 7, 8, 5}; Plane Surface(5) = {5}; Surface Loop(1) = {3, 2, 1, 4, 5}; Volume(1) = {1}; Line(95) = {17, 18}; Line(96) = {18, 19}; Line(97) = {19, 20}; Line(98) = {20, 17}; Line(100) = {13, 14}; Line(101) = {13, 17}; Line(102) = {18, 14}; Line(103) = {14, 15}; Line(104) = {19, 15}; Line(105) = {15, 16}; Line(106) = {20, 16}; Line(107) = {16, 13}; Line Loop(7) = {95, 96, 97, 98}; Plane Surface(8) = {7}; Line Loop(8) = {100, 103, 105, 107}; Plane Surface(9) = {8}; Line Loop(9) = {-106, 98, -101, -107}; Plane Surface(10) = {9}; Line Loop(10) = {101, 95, 102, -100}; Plane Surface(11) = {10}; Line Loop(11) = {-102, 96, 104, -103}; Plane Surface(12) = {11}; Line Loop(12) = {-104, 97, 106, -105}; Plane Surface(13) = {12}; Surface Loop(101) = {8, 9, 10, 11, 12, 13}; Volume(3) = {101}; Physical Surface("top_surface") = {5}; -Physical Surface("flat") = {9}; +Physical Surface("contact_bottom") = {9}; Physical Surface("bottom_surface") = {8}; Physical Volume("top_body") = {1}; Physical Volume("bot_body") = {3}; -Physical Surface("curved") = {1, 2, 4, 3}; +Physical Surface("contact_top") = {1, 2, 4, 3}; diff --git a/test/test_model/test_contact_mechanics_model/material.dat b/test/test_model/test_contact_mechanics_model/material.dat index b641cb2cb..4114e603f 100644 --- a/test/test_model/test_contact_mechanics_model/material.dat +++ b/test/test_model/test_contact_mechanics_model/material.dat @@ -1,25 +1,26 @@ material elastic [ name = bot_body rho = 7800 # density - E = 210 # young's modulu - nu = 0.3 # poisson's ratio + E = 1e5 # young's modulu + nu = 0. # poisson's ratio ] material elastic [ name = top_body rho = 7800 # density - E = 2.1e11 # young's modulu - nu = 0.3 # poisson's ratio + E = 1 # young's modulu + nu = 0. # poisson's ratio ] contact_detector [ type = implicit - master = flat - slave = curved + master = contact_bottom + slave = contact_top ] contact_resolution penalty [ - name = curved - epsilon = 100 + name = contact_top + mu = 0.0 + epsilon = 1e6 epsilon_t = 100 -] \ No newline at end of file +] \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/test_circles_2d.cc b/test/test_model/test_contact_mechanics_model/test_circles_2d.cc new file mode 100644 index 000000000..10b33cb5f --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_circles_2d.cc @@ -0,0 +1,98 @@ +/** + * @file test_contact_mechanics_model.cc + * + * @author Mohit Pundir + * + * @date creation: Tue Apr 30 2019 + * @date last modification: Tue Apr 30 2019 + * + * @brief Test for contact mechanics model 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 "solid_mechanics_model.hh" +#include "contact_mechanics_model.hh" +#include "coupler_solid_contact.hh" +#include "non_linear_solver.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + + const UInt spatial_dimension = 2; + initialize("material.dat", argc, argv); + + auto increment = 1e-3; + auto nb_steps = 10; + + Mesh mesh(spatial_dimension); + mesh.read("circles_2d.msh"); + + CouplerSolidContact coupler(mesh); + + auto & solid = coupler.getSolidMechanicsModel(); + auto & contact = coupler.getContactMechanicsModel(); + + solid.initFull( _analysis_method = _static); + contact.initFull(_analysis_method = _implicit_contact); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bot_body"); + //solid.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "bot_body"); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "top"); + + coupler.initFull(_analysis_method = _implicit_contact); + + auto & solver = coupler.getNonLinearSolver(); + solver.set("max_iterations", 50); + solver.set("threshold", 1e-8); + solver.set("convergence_type", _scc_residual); + + coupler.setBaseName("test-circles-2d"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("normals"); + coupler.addDumpFieldVector("contact_force"); + coupler.addDumpFieldVector("external_force"); + coupler.addDumpFieldVector("internal_force"); + coupler.addDumpField("gaps"); + coupler.addDumpField("blocked_dofs"); + coupler.addDumpField("grad_u"); + coupler.addDumpField("stress"); + + for (auto i : arange(nb_steps)) { + + solid.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "bot_body"); + coupler.solveStep(); + + coupler.dump(); + + std::cerr << "Step " << i << std::endl; + } + + finalize(); + return EXIT_SUCCESS; +} + diff --git a/test/test_model/test_contact_mechanics_model/test_contact_mechanics_model.cc b/test/test_model/test_contact_mechanics_model/test_contact_mechanics_model.cc deleted file mode 100644 index 79d2e21f2..000000000 --- a/test/test_model/test_contact_mechanics_model/test_contact_mechanics_model.cc +++ /dev/null @@ -1,60 +0,0 @@ -/** - * @file test_contact_mechanics_model.cc - * - * @author Mohit Pundir - * - * @date creation: Tue Apr 30 2019 - * @date last modification: Tue Apr 30 2019 - * - * @brief Test for contact mechanics model 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 "contact_mechanics_model.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -int main(int argc, char *argv[]) { - - initialize("material.dat", argc, argv); - const UInt spatial_dimension = 3; - - Mesh mesh(spatial_dimension); - mesh.read("contact_3d.msh"); - - ContactMechanicsModel model(mesh); - model.initFull(_analysis_method = _implicit_contact); - - model.setBaseNameToDumper("contact_mechanics", "contact"); - model.addDumpFieldVectorToDumper("contact_mechanics", "contact_force"); - model.addDumpFieldVectorToDumper("contact_mechanics", "external_force"); - model.addDumpFieldVectorToDumper("contact_mechanics", "normals"); - model.addDumpFieldToDumper("contact_mechanics", "gaps"); - model.addDumpFieldToDumper("contact_mechanics", "areas"); - - model.search(); - model.dump("contact_mechanics"); - - return 0; -} - diff --git a/test/test_model/test_contact_mechanics_model/test_contact_mechanics_model_2d.cc b/test/test_model/test_contact_mechanics_model/test_contact_mechanics_model_2d.cc deleted file mode 100644 index 3f18c85a8..000000000 --- a/test/test_model/test_contact_mechanics_model/test_contact_mechanics_model_2d.cc +++ /dev/null @@ -1,59 +0,0 @@ -/** - * @file test_contact_mechanics_model.cc - * - * @author Mohit Pundir - * - * @date creation: Tue Apr 30 2019 - * @date last modification: Tue Apr 30 2019 - * - * @brief Test for contact mechanics model 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 "contact_mechanics_model.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -int main(int argc, char *argv[]) { - - initialize("material.dat", argc, argv); - const UInt spatial_dimension = 2; - - Mesh mesh(spatial_dimension); - mesh.read("contact_2d.msh"); - - ContactMechanicsModel model(mesh); - model.initFull(_analysis_method = _implicit_contact); - - model.setBaseNameToDumper("contact_mechanics", "contact-2d"); - model.addDumpFieldVectorToDumper("contact_mechanics", "contact_force"); - model.addDumpFieldVectorToDumper("contact_mechanics", "normals"); - model.addDumpFieldToDumper("contact_mechanics", "gaps"); - model.addDumpFieldToDumper("contact_mechanics", "areas"); - - model.search(); - model.dump("contact_mechanics"); - - return 0; -} - diff --git a/test/test_model/test_contact_mechanics_model/test_coupler/CMakeLists.txt b/test/test_model/test_contact_mechanics_model/test_coupler/CMakeLists.txt index a57670c2e..45f11aae3 100644 --- a/test/test_model/test_contact_mechanics_model/test_coupler/CMakeLists.txt +++ b/test/test_model/test_contact_mechanics_model/test_coupler/CMakeLists.txt @@ -1,32 +1,42 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Tue May 07 2019 # @date last modification: Tue May 07 2019 # # @brief configuration for testing coupler class for 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 . # # @section DESCRIPTION # #=============================================================================== add_mesh (coupling coupling.geo 2 1) +add_mesh (flat_on_flat flat_on_flat.geo 2 1) + register_test( test_contact_coupling SOURCES test_contact_coupling.cc DEPENDS coupling FILES_TO_COPY material.dat PACKAGE contact_mechanics model_couplers ) + + +register_test( test_coupled_stiffness + SOURCES test_coupled_stiffness.cc + DEPENDS flat_on_flat + FILES_TO_COPY material_stiffness.dat + PACKAGE contact_mechanics model_couplers + ) diff --git a/test/test_model/test_contact_mechanics_model/test_coupler/flat_on_flat.geo b/test/test_model/test_contact_mechanics_model/test_coupler/flat_on_flat.geo new file mode 100644 index 000000000..5085f88dd --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_coupler/flat_on_flat.geo @@ -0,0 +1,33 @@ +a = 1.0; +r = a; +gap = 0; +nbpoints = 1.0; +cl1 = a/(nbpoints); +cl2 = r/(nbpoints); + +Point(1) = {-a/2, 0 + gap, 0, cl1}; +Point(2) = {a/2, 0 + gap, 0, cl1}; +Point(3) = {a/2, a/2., 0, cl1}; +Point(4) = {-a/2, a/2., 0, cl1}; +Point(5) = {r/2, 0, 0, cl2}; +Point(6) = {-r/2, 0, 0, cl2}; +Point(7) = {-r/2, -r/2, 0, 10*cl2}; +Point(8) = { r/2, -r/2, 0, 10*cl2}; +Line(2) = {1, 2}; +Line(3) = {2, 3}; +Line(4) = {3, 4}; +Line(5) = {4, 1}; +Line(6) = {7, 8}; +Line(7) = {8, 5}; +Line(8) = {5, 6}; +Line(9) = {6, 7}; +Line Loop(1) = {2, 3, 4, 5}; +Line Loop(2) = {6, 7, 8, 9}; +Plane Surface(1) = {1}; +Plane Surface(2) = {2}; +Physical Line("contact_top") = {2}; +Physical Line("contact_bottom") = {8}; +Physical Line("top") = {4}; +Physical Line("bottom") = {6}; +Physical Surface("bot_body") = {2}; +Physical Surface("top_body") = {1}; \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/material.dat b/test/test_model/test_contact_mechanics_model/test_coupler/material_stiffness.dat similarity index 66% copy from test/test_model/test_contact_mechanics_model/material.dat copy to test/test_model/test_contact_mechanics_model/test_coupler/material_stiffness.dat index b641cb2cb..68376f8e8 100644 --- a/test/test_model/test_contact_mechanics_model/material.dat +++ b/test/test_model/test_contact_mechanics_model/test_coupler/material_stiffness.dat @@ -1,25 +1,24 @@ material elastic [ name = bot_body rho = 7800 # density - E = 210 # young's modulu + E = 1. # young's modulu nu = 0.3 # poisson's ratio ] material elastic [ name = top_body rho = 7800 # density - E = 2.1e11 # young's modulu + E = 1e5 # young's modulu nu = 0.3 # poisson's ratio ] contact_detector [ type = implicit - master = flat - slave = curved + master = contact_top + slave = contact_bottom ] contact_resolution penalty [ - name = curved - epsilon = 100 - epsilon_t = 100 + name = contact_bottom + epsilon = 1e6 ] \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/test_coupler/test_coupled_stiffness.cc b/test/test_model/test_contact_mechanics_model/test_coupler/test_coupled_stiffness.cc new file mode 100644 index 000000000..d8ccc88fb --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_coupler/test_coupled_stiffness.cc @@ -0,0 +1,130 @@ +/** + * @file test_contact_mechanics_model.cc + * + * @author Mohit Pundir + * + * @date creation: Tue Apr 30 2019 + * @date last modification: Tue Apr 30 2019 + * + * @brief Test for contact mechanics model 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 "solid_mechanics_model.hh" +#include "contact_mechanics_model.hh" +#include "coupler_solid_contact.hh" +#include "non_linear_solver.hh" +#include "sparse_matrix.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + + const UInt spatial_dimension = 2; + initialize("material_stiffness.dat", argc, argv); + + Mesh mesh(spatial_dimension); + mesh.read("flat_on_flat.msh"); + + CouplerSolidContact coupler(mesh); + + auto & solid = coupler.getSolidMechanicsModel(); + auto & contact = coupler.getContactMechanicsModel(); + + solid.initFull( _analysis_method = _static); + contact.initFull(_analysis_method = _implicit_contact); + + 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"); + + coupler.initFull(_analysis_method = _implicit_contact); + + auto & solver = coupler.getNonLinearSolver(); + solver.set("max_iterations", 1); + solver.set("threshold", 1e-1); + solver.set("convergence_type", _scc_solution); + + coupler.setBaseName("test-coupled-stiffness"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("normals"); + coupler.addDumpFieldVector("contact_force"); + coupler.addDumpFieldVector("external_force"); + coupler.addDumpFieldVector("internal_force"); + coupler.addDumpField("gaps"); + coupler.addDumpField("blocked_dofs"); + coupler.addDumpField("grad_u"); + coupler.addDumpField("stress"); + + + auto & before_assembly = + const_cast(coupler.getDOFManager().getNewMatrix("K", _symmetric)); + + solid.assembleStiffnessMatrix(); + + auto & solid_assembly = + const_cast(coupler.getDOFManager().getMatrix("K")); + solid_assembly.saveMatrix("solid_assembly.mtx"); + + auto & displacement = solid.getDisplacement(); + + contact.search(displacement); + contact.assembleStiffnessMatrix(); + + auto contact_map = contact.getContactMap(); + auto nb_contacts = contact_map.size(); + + auto & contact_assembly = + const_cast(coupler.getDOFManager().getMatrix("K")); + contact_assembly.saveMatrix("contact_assembly.mtx"); + + solid.assembleInternalForces(); + contact.assembleInternalForces(); + + coupler.dump(); + + Array & contact_force = contact.getInternalForce(); + + for (UInt n : arange(contact_force.size())) { + std::cerr << contact_force(n, 1) << std::endl; + } + + if (solid_assembly.size() == contact_assembly.size() and nb_contacts > 0) { + std::cerr << "size of stiffness matrix of solid = " << solid_assembly.size() << std::endl; + std::cerr << "size of stiffness matrix of coupled = " << contact_assembly.size() << std::endl; + std::cerr << "number of contacts = " << nb_contacts << std::endl; + + for (auto & pair : contact_map) { + std::cerr << "Node " << pair.first << " in contact with " << pair.second.master << + " of gap " << pair.second.gap << std::endl; + } + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} + diff --git a/test/test_model/test_contact_mechanics_model/test_detector/CMakeLists.txt b/test/test_model/test_contact_mechanics_model/test_detector/CMakeLists.txt index 8c147f00d..dcf2508b6 100644 --- a/test/test_model/test_contact_mechanics_model/test_detector/CMakeLists.txt +++ b/test/test_model/test_contact_mechanics_model/test_detector/CMakeLists.txt @@ -1,50 +1,59 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Wed Dec 18 2018 # @date last modification: Wed Apr 29 2019 # # @brief configuration for contact detection tests # # @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 . # # @section DESCRIPTION # #=============================================================================== add_mesh (implicit_2d implicit_2d.geo 2 1) add_mesh (spheres_2d spheres_2d.geo 2 1) add_mesh (implicit_3d implicit_3d.geo 3 1) - +add_mesh (paraview_2d paraview_2d.geo 2 1) register_test(test_implicit_3d SOURCES test_detection_implicit_3d.cc DEPENDS implicit_3d FILES_TO_COPY options.dat PACKAGE contact_mechanics ) register_test(test_implicit_2d SOURCES test_detection_implicit_2d.cc DEPENDS implicit_2d FILES_TO_COPY options.dat PACKAGE contact_mechanics ) + +register_test(test_paraview + SOURCES test_paraview.cc + DEPENDS paraview_2d + FILES_TO_COPY material.dat + PACKAGE contact_mechanics + ) + + register_test(test_spheres_2d SOURCES test_detection_implicit_2d.cc DEPENDS spheres_2d FILES_TO_COPY options.dat PACKAGE contact_mechanics ) diff --git a/test/test_model/test_contact_mechanics_model/material.dat b/test/test_model/test_contact_mechanics_model/test_detector/material.dat similarity index 50% copy from test/test_model/test_contact_mechanics_model/material.dat copy to test/test_model/test_contact_mechanics_model/test_detector/material.dat index b641cb2cb..4114e603f 100644 --- a/test/test_model/test_contact_mechanics_model/material.dat +++ b/test/test_model/test_contact_mechanics_model/test_detector/material.dat @@ -1,25 +1,26 @@ material elastic [ name = bot_body rho = 7800 # density - E = 210 # young's modulu - nu = 0.3 # poisson's ratio + E = 1e5 # young's modulu + nu = 0. # poisson's ratio ] material elastic [ name = top_body rho = 7800 # density - E = 2.1e11 # young's modulu - nu = 0.3 # poisson's ratio + E = 1 # young's modulu + nu = 0. # poisson's ratio ] contact_detector [ type = implicit - master = flat - slave = curved + master = contact_bottom + slave = contact_top ] contact_resolution penalty [ - name = curved - epsilon = 100 + name = contact_top + mu = 0.0 + epsilon = 1e6 epsilon_t = 100 -] \ No newline at end of file +] \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/circles_2d.geo b/test/test_model/test_contact_mechanics_model/test_detector/paraview_2d.geo similarity index 65% copy from test/test_model/test_contact_mechanics_model/circles_2d.geo copy to test/test_model/test_contact_mechanics_model/test_detector/paraview_2d.geo index f29094576..ec7b08464 100644 --- a/test/test_model/test_contact_mechanics_model/circles_2d.geo +++ b/test/test_model/test_contact_mechanics_model/test_detector/paraview_2d.geo @@ -1,43 +1,43 @@ -cl1 = 0.0025; -cl2 = 0.0025; +cl1 = 0.001; +cl2 = 0.005; cl3 = 0.005; Dy = 0.0; radius = 0.1; y = 0.1; -epsilon = 0.0; -Point(1) = {0, y, 0, cl1}; -Point(2) = {radius, radius + y, 0, cl2}; -Point(3) = {-radius, radius + y, 0, cl2}; +epsilon = 1e-10; +Point(1) = {0, y + epsilon, 0, cl1}; +Point(2) = {radius, radius + y + epsilon, 0, cl2}; +Point(3) = {-radius, radius + y + epsilon, 0, cl2}; -Point(11) = {0, y, 0, cl1}; +Point(11) = {0, y, 0, cl2}; Point(12) = {radius, -radius + y, 0, cl2}; Point(13) = {-radius,- radius + y, 0, cl2}; Point(8) = {0, radius + y, 0, cl2}; Point(18) = {0, -radius + y, 0, cl2}; Circle(1) = {3, 8, 1}; Circle(2) = {1, 8, 2}; Circle(11) = {13, 18, 11}; Circle(12) = {11, 18, 12}; Line(3) = {2, 8}; Line(13) = {8, 3}; Line(30) = {12, 18}; Line(31) = {18, 13}; Line Loop(9) = {2, 3, 13, 1}; Plane Surface(9) = {9}; Line Loop(19) = {-12, -30, -31, -11}; Plane Surface(19) = {19}; -Physical Line("rigid") = {11, 12}; -Physical Line("contact_surface") = {1, 2}; +Physical Line("contact_bottom") = {11, 12}; +Physical Line("contact_top") = {1, 2}; Physical Line("bottom") = {30, 31}; Physical Line("top") = {3, 13}; Physical Surface("bot_body") = {19}; Physical Surface("top_body") = {9}; \ No newline at end of file diff --git a/test/test_model/test_contact_mechanics_model/test_detector/test_paraview.cc b/test/test_model/test_contact_mechanics_model/test_detector/test_paraview.cc new file mode 100644 index 000000000..abfa4955e --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_detector/test_paraview.cc @@ -0,0 +1,111 @@ +/** + * @file test_contact_mechanics_model.cc + * + * @author Mohit Pundir + * + * @date creation: Tue Apr 30 2019 + * @date last modification: Tue Apr 30 2019 + * + * @brief Test for contact mechanics model 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 "solid_mechanics_model.hh" +#include "contact_mechanics_model.hh" +#include "coupler_solid_contact.hh" +#include "non_linear_solver.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + + const UInt spatial_dimension = 2; + initialize("material.dat", argc, argv); + + auto increment = 1e-3; + auto nb_steps = 10; + + Mesh mesh(spatial_dimension); + mesh.read("paraview_2d.msh"); + + CouplerSolidContact coupler(mesh); + + auto & solid = coupler.getSolidMechanicsModel(); + auto & contact = coupler.getContactMechanicsModel(); + + solid.initFull( _analysis_method = _static); + contact.initFull(_analysis_method = _implicit_contact); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bot_body"); + solid.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "bot_body"); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "top"); + + coupler.initFull(_analysis_method = _implicit_contact); + + coupler.setBaseName("test-circles-2d"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("normals"); + coupler.addDumpFieldVector("contact_force"); + coupler.addDumpFieldVector("external_force"); + coupler.addDumpFieldVector("internal_force"); + coupler.addDumpField("gaps"); + coupler.addDumpField("blocked_dofs"); + coupler.addDumpField("grad_u"); + coupler.addDumpField("stress"); + + auto & before_assembly = + const_cast(coupler.getDOFManager().getNewMatrix("K", _symmetric)); + + solid.assembleStiffnessMatrix(); + + auto & solid_assembly = + const_cast(coupler.getDOFManager().getMatrix("K")); + solid_assembly.saveMatrix("solid_assembly.mtx"); + + auto & displacement = solid.getDisplacement(); + + contact.search(displacement); + contact.assembleStiffnessMatrix(); + + auto contact_map = contact.getContactMap(); + auto nb_contacts = contact_map.size(); + + auto & contact_assembly = + const_cast(coupler.getDOFManager().getMatrix("K")); + contact_assembly.saveMatrix("contact_assembly.mtx"); + + solid.assembleInternalForces(); + contact.assembleInternalForces(); + + coupler.dump(); + + + + finalize(); + return EXIT_SUCCESS; +} + diff --git a/test/test_model/test_contact_mechanics_model/test_flat_2d.cc b/test/test_model/test_contact_mechanics_model/test_flat_2d.cc new file mode 100644 index 000000000..07bee6731 --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_flat_2d.cc @@ -0,0 +1,91 @@ +/** + * @file test_flat_2d.cc + * + * @author Mohit Pundir + * + * @date creation: Tue Apr 30 2019 + * @date last modification: Tue Apr 30 2019 + * + * @brief Test for contact mechanics model 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 "solid_mechanics_model.hh" +#include "contact_mechanics_model.hh" +#include "coupler_solid_contact.hh" +#include "non_linear_solver.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + + const UInt spatial_dimension = 2; + initialize("material.dat", argc, argv); + + Mesh mesh(spatial_dimension); + mesh.read("flat_2d.msh"); + + CouplerSolidContact coupler(mesh); + + auto & solid = coupler.getSolidMechanicsModel(); + auto & contact = coupler.getContactMechanicsModel(); + + solid.initFull( _analysis_method = _static); + contact.initFull(_analysis_method = _implicit_contact); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bot_body"); + solid.applyBC(BC::Dirichlet::IncrementValue(0.005, _y), "bot_body"); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "top"); + + coupler.initFull(_analysis_method = _implicit_contact); + + auto & solver = coupler.getNonLinearSolver(); + solver.set("max_iterations", 1000); + solver.set("threshold", 1e-8); + solver.set("convergence_type", _scc_solution); + + coupler.setBaseName("test-flat-2d"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("normals"); + coupler.addDumpFieldVector("contact_force"); + coupler.addDumpFieldVector("external_force"); + coupler.addDumpFieldVector("internal_force"); + coupler.addDumpField("gaps"); + coupler.addDumpField("blocked_dofs"); + coupler.addDumpField("grad_u"); + coupler.addDumpField("stress"); + + coupler.dump(); + + coupler.solveStep(); + + coupler.dump(); + + finalize(); + return EXIT_SUCCESS; +} + diff --git a/test/test_model/test_contact_mechanics_model/test_hertz_2d.cc b/test/test_model/test_contact_mechanics_model/test_hertz_2d.cc new file mode 100644 index 000000000..23145134f --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_hertz_2d.cc @@ -0,0 +1,99 @@ +/** + * @file test_contact_mechanics_model.cc + * + * @author Mohit Pundir + * + * @date creation: Tue Apr 30 2019 + * @date last modification: Tue Apr 30 2019 + * + * @brief Test for contact mechanics model 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 "solid_mechanics_model.hh" +#include "contact_mechanics_model.hh" +#include "coupler_solid_contact.hh" +#include "non_linear_solver.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + + const UInt spatial_dimension = 2; + initialize("material.dat", argc, argv); + + auto increment = 1e-3; + auto nb_steps = 1; + + Mesh mesh(spatial_dimension); + mesh.read("hertz_2d.msh"); + + CouplerSolidContact coupler(mesh); + + auto & solid = coupler.getSolidMechanicsModel(); + auto & contact = coupler.getContactMechanicsModel(); + + solid.initFull( _analysis_method = _static); + contact.initFull(_analysis_method = _implicit_contact); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bot_body"); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "top"); + + coupler.initFull(_analysis_method = _implicit_contact); + + auto & solver = coupler.getNonLinearSolver(); + solver.set("max_iterations", 1000); + solver.set("threshold", 1e-8); + solver.set("convergence_type", _scc_residual); + + coupler.setBaseName("trial-hertz-2d"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("normals"); + coupler.addDumpFieldVector("contact_force"); + coupler.addDumpFieldVector("external_force"); + coupler.addDumpFieldVector("internal_force"); + coupler.addDumpField("gaps"); + coupler.addDumpField("blocked_dofs"); + coupler.addDumpField("grad_u"); + coupler.addDumpField("stress"); + + coupler.dump(); + + for (auto i : arange(nb_steps)) { + + solid.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "bot_body"); + coupler.solveStep(); + + coupler.dump(); + + std::cerr << "Step " << i << std::endl; + } + + finalize(); + return EXIT_SUCCESS; +} + diff --git a/test/test_model/test_contact_mechanics_model/test_hertz_3d.cc b/test/test_model/test_contact_mechanics_model/test_hertz_3d.cc new file mode 100644 index 000000000..346cb3a9d --- /dev/null +++ b/test/test_model/test_contact_mechanics_model/test_hertz_3d.cc @@ -0,0 +1,100 @@ +/** + * @file test_contact_mechanics_model.cc + * + * @author Mohit Pundir + * + * @date creation: Tue Apr 30 2019 + * @date last modification: Tue Apr 30 2019 + * + * @brief Test for contact mechanics model 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 "solid_mechanics_model.hh" +#include "contact_mechanics_model.hh" +#include "coupler_solid_contact.hh" +#include "non_linear_solver.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + + const UInt spatial_dimension = 3; + initialize("material.dat", argc, argv); + + auto increment = 5e-3; + auto nb_steps = 10; + + Mesh mesh(spatial_dimension); + mesh.read("hertz_3d.msh"); + + CouplerSolidContact coupler(mesh); + + auto & solid = coupler.getSolidMechanicsModel(); + auto & contact = coupler.getContactMechanicsModel(); + + solid.initFull( _analysis_method = _static); + contact.initFull(_analysis_method = _implicit_contact); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bot_body"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "bot_body"); + + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top_surface"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "top_surface"); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "top_surface"); + + coupler.initFull(_analysis_method = _implicit_contact); + + auto & solver = coupler.getNonLinearSolver(); + solver.set("max_iterations", 1000); + solver.set("threshold", 1e-8); + solver.set("convergence_type", _scc_residual); + + coupler.setBaseName("test-hertz-3d"); + coupler.addDumpFieldVector("displacement"); + coupler.addDumpFieldVector("normals"); + coupler.addDumpFieldVector("contact_force"); + coupler.addDumpFieldVector("external_force"); + coupler.addDumpFieldVector("internal_force"); + coupler.addDumpField("gaps"); + coupler.addDumpField("blocked_dofs"); + coupler.addDumpField("grad_u"); + coupler.addDumpField("stress"); + + coupler.dump(); + + for (auto i : arange(nb_steps)) { + + solid.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "bot_body"); + coupler.solveStep(); + coupler.dump(); + + std::cerr << "Step " << i << std::endl; + } + + finalize(); + return EXIT_SUCCESS; +} +