diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index 7fc7ef6dd..27a888b70 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,614 +1,282 @@ /** * @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, const MemoryID & memory_id) : ContactDetector(mesh, mesh.getNodes(), id, memory_id) {} /* -------------------------------------------------------------------------- */ ContactDetector::ContactDetector(Mesh & mesh, Array positions, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), Parsable(ParserType::_contact_detector, id), mesh(mesh), positions(0, mesh.getSpatialDimension()) { AKANTU_DEBUG_IN(); this->spatial_dimension = mesh.getSpatialDimension(); this->positions.copy(positions); this->parseSection(); 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); } } /* -------------------------------------------------------------------------- */ void ContactDetector::search(Array & elements, Array & gaps, Array & normals, Array & projections) { UInt surface_dimension = spatial_dimension - 1; this->mesh.fillNodesToElements(surface_dimension); this->computeMaximalDetectionDistance(); contact_pairs.clear(); SpatialGrid master_grid(spatial_dimension); SpatialGrid slave_grid(spatial_dimension); this->globalSearch(slave_grid, master_grid); this->localSearch(slave_grid, master_grid); this->createContactElements(elements, gaps, normals, projections); } /* -------------------------------------------------------------------------- */ void ContactDetector::globalSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid) { auto & master_list = surface_selector->getMasterList(); auto & slave_list = surface_selector->getSlaveList(); 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 /// 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)) { /// check for self contact if (slave_node == master_node) continue; bool is_valid = true; Array elements; this->mesh.getAssociatedElements(slave_node, elements); for (auto & elem : elements) { if (elem.kind() != _ek_regular) continue; Vector connectivity = const_cast(this->mesh).getConnectivity(elem); auto node_iter = std::find(connectivity.begin(), connectivity.end(), master_node); if (node_iter != connectivity.end()) { is_valid = false; break; } } 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 and is_valid) { 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; - - std::map previous_contact_map = contact_map; - contact_map.clear(); - - auto get_connectivity = [&](auto & slave, auto & master) { - Vector master_conn = - const_cast(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; - }; - - for (auto & pairs : contact_pairs) { - - const auto & slave_node = pairs.first; - const auto & master_node = pairs.second; - - Array all_elements; - this->mesh.getAssociatedElements(master_node, all_elements); - - Array boundary_elements; - this->filterBoundaryElements(all_elements, boundary_elements); - - Array gaps(boundary_elements.size(), 1, "gaps"); - Array normals(boundary_elements.size(), spatial_dimension, "normals"); - Array projections(boundary_elements.size(), surface_dimension, - "projections"); - - auto index = this->computeOrthogonalProjection( - slave_node, boundary_elements, normals, gaps, projections); - - if (index == UInt(-1)) { - continue; - } - - auto connectivity = get_connectivity(slave_node, boundary_elements[index]); - - // assign contact element attributes - contact_map[slave_node].setMaster(boundary_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(surface_dimension)[index], true)); - contact_map[slave_node].setConnectivity(connectivity); - - // tangent computation on master surface - Matrix tangents(surface_dimension, spatial_dimension); - this->computeTangentsOnElement(contact_map[slave_node].master, - contact_map[slave_node].projection, - contact_map[slave_node].normal, tangents); - - contact_map[slave_node].setTangent(tangents); - - // friction calculation requires history of previous natural - // projection as well as traction - auto search = previous_contact_map.find(slave_node); - if (search != previous_contact_map.end()) { - auto previous_projection = - previous_contact_map[slave_node].getPreviousProjection(); - contact_map[slave_node].setPreviousProjection(previous_projection); - - } else { - Vector previous_projection(surface_dimension, 0.); - contact_map[slave_node].setPreviousProjection(previous_projection); - } - - // to ensure the self contact between surface does not lead to - // detection of master element which is closet but not - // orthogonally opposite to the slave surface - bool is_valid_self_contact = - this->checkValidityOfSelfContact(slave_node, contact_map[slave_node]); - - if (!is_valid_self_contact) { - contact_map.erase(slave_node); - } - } - - contact_pairs.clear(); - }*/ - /* -------------------------------------------------------------------------- */ void ContactDetector::createContactElements(Array & contact_elements, Array & gaps, Array & normals, Array & projections) { auto surface_dimension = spatial_dimension - 1; Real alpha; switch (detection_type) { case _explicit: { alpha = 1.0; break; } case _implicit: { alpha = -1.0; break; } default: AKANTU_EXCEPTION(detection_type << " is not a valid contact detection type"); break; } for (auto & pairs : contact_pairs) { const auto & slave_node = pairs.first; Vector slave(spatial_dimension); for (UInt s : arange(spatial_dimension)) slave(s) = this->positions(slave_node, s); const auto & master_node = pairs.second; Array elements; this->mesh.getAssociatedElements(master_node, elements); auto & gap = gaps.begin()[slave_node]; Vector normal(normals.begin(spatial_dimension)[slave_node]); Vector projection(projections.begin(surface_dimension)[slave_node]); auto index = GeometryUtils::orthogonalProjection(mesh, positions, slave, elements, gap, projection, normal, alpha); // if not a valid projection is found on patch of elements index is -1 if (index == UInt(-1)) continue; // if not a valid self contact if (!isValidSelfContact(slave_node, gap, normal)) continue; // create contact element contact_elements.push_back(ContactElement(slave_node, elements[index])); } contact_pairs.clear(); } - -/* -------------------------------------------------------------------------- */ -/*UInt 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); - - UInt counter = 0; - UInt index = UInt(-1); - - Real min_gap = std::numeric_limits::max(); - - 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); - - gap = this->computeGap(query, real_projection); - - // check if gap is valid or not - // to check this we need normal on master element, vector from - // real projection to slave node, if it is explicit detection - // scheme, we want it to interpenetrate, the dot product should be - // negative and -1.0 and for implciit detection , the dot product - // should be positive and 1.0 - - bool is_valid = this->checkValidityOfProjection(projection); - - auto master_to_slave = query - real_projection; - auto norm = master_to_slave.norm(); - if (norm != 0) - master_to_slave /= norm; - - Real tolerance = 1e-8; - - switch (detection_type) { - case _explicit: { - auto product = master_to_slave.dot(normal); - auto variation = std::abs(product + 1.0); - - if (variation <= tolerance and gap <= min_gap and is_valid) { - min_gap = gap; - index = counter; - gap *= -1.0; - } - break; - } - case _implicit: { - auto product = master_to_slave.dot(normal); - auto variation = std::abs(product - 1.0); - - if (variation <= tolerance and gap <= min_gap and is_valid) { - min_gap = gap; - index = counter; - gap *= -1.0; - } - } - default: - break; - } - - counter++; - } - - return index; - }*/ - -/* -------------------------------------------------------------------------- */ -/*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::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::computeTangentsOnElement(const Element & el, - Vector & projection, - Vector & normal, - 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(); - - // to ensure that direction of tangents are correct, cross product - // of tangents should give the normal vector computed earlier - switch (spatial_dimension) { - case 2: { - Vector e_z(3); - e_z[0] = 0.; - e_z[1] = 0.; - e_z[2] = 1.; - - Vector tangent(3); - tangent[0] = tangents(0, 0); - tangent[1] = tangents(0, 1); - tangent[2] = 0.; - - auto exp_normal = e_z.crossProduct(tangent); - - auto & cal_normal = normal; - - auto ddot = cal_normal.dot(exp_normal); - if (ddot < 0) { - tangents *= -1.0; - } - - break; - } - case 3: { - auto tang_trans = tangents.transpose(); - auto tang1 = Vector(tang_trans(0)); - auto tang2 = Vector(tang_trans(1)); - - auto tang1_cross_tang2 = tang1.crossProduct(tang2); - auto exp_normal = tang1_cross_tang2 / tang1_cross_tang2.norm(); - - auto & cal_normal = normal; - - auto ddot = cal_normal.dot(exp_normal); - if (ddot < 0) { - tang_trans(1) *= -1.0; - } - - tangents = tang_trans.transpose(); - break; - } - default: - break; - } - }*/ - } // namespace akantu diff --git a/src/model/contact_mechanics/contact_detector_inline_impl.cc b/src/model/contact_mechanics/contact_detector_inline_impl.cc index 81cc68b7c..4bd22b2e2 100644 --- a/src/model/contact_mechanics/contact_detector_inline_impl.cc +++ b/src/model/contact_mechanics/contact_detector_inline_impl.cc @@ -1,403 +1,346 @@ /** * @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 tolerance = 1e-3; for (auto xi : projection) { if (xi >= -1.0 - tolerance and xi <= 1.0 + tolerance) 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; } 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); } /* -------------------------------------------------------------------------- */ inline void ContactDetector::computeMaximalDetectionDistance() { AKANTU_DEBUG_IN(); Real elem_size; Real max_elem_size = std::numeric_limits::min(); Real min_elem_size = std::numeric_limits::max(); auto & master_nodes = this->surface_selector->getMasterList(); for (auto & master : master_nodes) { Array elements; this->mesh.getAssociatedElements(master, elements); for (auto element : elements) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element.type); Matrix elem_coords(spatial_dimension, nb_nodes_per_element); this->coordinatesOfElement(element, elem_coords); elem_size = FEEngine::getElementInradius(elem_coords, element.type); max_elem_size = std::max(max_elem_size, elem_size); min_elem_size = std::min(min_elem_size, elem_size); } } AKANTU_DEBUG_INFO("The maximum element size : " << max_elem_size); this->min_dd = min_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 = const_cast(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); } } // to ensure that normal is always outwards from master surface const auto & element_to_subelement = mesh.getElementToSubelement(element.type)(element.element); Vector outside(spatial_dimension); mesh.getBarycenter(element, outside); // check if mesh facets exists for cohesive elements contact Vector inside(spatial_dimension); if (mesh.isMeshFacets()) { mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); } else { mesh.getBarycenter(element_to_subelement[0], inside); } Vector inside_to_outside = outside - inside; auto projection = inside_to_outside.dot(normal); if (projection < 0) { normal *= -1.0; } } /* -------------------------------------------------------------------------- */ 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); } } } /* -------------------------------------------------------------------------- */ inline Real ContactDetector::computeGap(Vector & slave, Vector & master) { Vector slave_to_master(spatial_dimension); slave_to_master = master - slave; Real gap = slave_to_master.norm(); return gap; } /* -------------------------------------------------------------------------- */ inline void ContactDetector::filterBoundaryElements(Array & elements, Array & boundary_elements) { for (auto elem : elements) { const auto & element_to_subelement = mesh.getElementToSubelement(elem.type)(elem.element); // for regular boundary elements if (element_to_subelement.size() == 1 and element_to_subelement[0].kind() == _ek_regular) { boundary_elements.push_back(elem); continue; } // for cohesive boundary elements UInt nb_subelements_regular = 0; for (auto subelem : element_to_subelement) { if (subelem == ElementNull) continue; if (subelem.kind() == _ek_regular) ++nb_subelements_regular; } auto nb_subelements = element_to_subelement.size(); if (nb_subelements_regular < nb_subelements) boundary_elements.push_back(elem); } } -/* -------------------------------------------------------------------------- */ -/*inline bool -ContactDetector::checkValidityOfSelfContact(const UInt & slave_node, - ContactElement & element) { - - UInt master_node; - - for (auto & pair : contact_pairs) { - if (pair.first == slave_node) { - master_node = pair.second; - break; - } - } - - Array elements; - this->mesh.getAssociatedElements(slave_node, elements); - - Vector slave_normal(spatial_dimension); - for (auto & elem : elements) { - if (elem.kind() != _ek_regular) - continue; - - Vector connectivity = - const_cast(this->mesh).getConnectivity(elem); - - Vector normal(spatial_dimension); - this->computeNormalOnElement(elem, normal); - - slave_normal = slave_normal + normal; - - auto node_iter = - std::find(connectivity.begin(), connectivity.end(), master_node); - if (node_iter != connectivity.end()) { - return false; - } - - - } - - if (std::abs(element.gap) > 2.0 * min_dd) { - return false; - } - - auto norm = slave_normal.norm(); - if (norm != 0) { - slave_normal /= norm; - } - - auto product = slave_normal.dot(element.normal); - if (product >= 0) { - return false; - } - - return true; - }*/ - - /* -------------------------------------------------------------------------- */ inline bool ContactDetector::isValidSelfContact(const UInt & slave_node, const Real & gap, const Vector & normal) { UInt master_node; // finding the master node corresponding to slave node for (auto & pair : contact_pairs) { if (pair.first == slave_node) { master_node = pair.second; break; } } Array slave_elements; this->mesh.getAssociatedElements(slave_node, slave_elements); // Check 1 : master node is not connected to elements connected to // slave node Vector slave_normal(spatial_dimension); for (auto & element : slave_elements) { if (element.kind() != _ek_regular) continue; Vector connectivity = const_cast(this->mesh).getConnectivity(element); // finding the normal at slave node by averaging of normals Vector normal(spatial_dimension); GeometryUtils::normal(mesh, positions, element, normal); slave_normal = slave_normal + normal; auto node_iter = std::find(connectivity.begin(), connectivity.end(), master_node); if (node_iter != connectivity.end()) return false; } // Check 2 : if gap is twice the size of smallest element if (std::abs(gap) > 2.0 * min_dd) return false; // Check 3 : check the directions of normal at slave node and at // master element, should be in opposite directions auto norm = slave_normal.norm(); if (norm != 0) slave_normal /= norm; auto product = slave_normal.dot(normal); if (product >= 0) return false; return true; } } // namespace 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 522c5d18a..ac9c6cd88 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,685 +1,698 @@ /** * @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" #include "group_manager_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id, memory_id) { 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, id + ":contact_detector"); + registerFEEngineObject("ContactFacetsFEEngine", mesh, + Model::spatial_dimension - 1); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::~ContactMechanicsModel() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // initalize the resolutions if (this->parser.getLastParsedFile() != "") { this->instantiateResolutions(); } this->initResolutions(); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::instantiateResolutions() { ParserSection model_section; bool is_empty; std::tie(model_section, is_empty) = this->getParserSection(); if (not is_empty) { auto model_resolutions = model_section.getSubSections(ParserType::_contact_resolution); for (const auto & section : model_resolutions) { this->registerNewResolution(section); } } auto sub_sections = this->parser.getSubSections(ParserType::_contact_resolution); for (const auto & section : sub_sections) { this->registerNewResolution(section); } if (resolutions.empty()) AKANTU_EXCEPTION("No contact resolutions where instantiated for the model" << getID()); are_resolutions_instantiated = true; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution(const ParserSection & section) { std::string res_name; std::string res_type = section.getName(); std::string opt_param = section.getOption(); try { std::string tmp = section.getParameter("name"); res_name = tmp; /** this can seem weird, but there is an ambiguous operator * overload that i couldn't solve. @todo remove the * weirdness of this code */ } catch (debug::Exception &) { AKANTU_ERROR("A contact resolution of type \'" << res_type << "\' in the input file has been defined without a name!"); } Resolution & res = this->registerNewResolution(res_name, res_type, opt_param); res.parseSection(section); return res; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution( const ID & res_name, const ID & res_type, const ID & opt_param) { AKANTU_DEBUG_ASSERT(resolutions_names_to_id.find(res_name) == resolutions_names_to_id.end(), "A resolution with this name '" << res_name << "' has already been registered. " << "Please use unique names for resolutions"); UInt res_count = resolutions.size(); resolutions_names_to_id[res_name] = res_count; std::stringstream sstr_res; sstr_res << this->id << ":" << res_count << ":" << res_type; ID res_id = sstr_res.str(); std::unique_ptr resolution = ResolutionFactory::getInstance().allocate(res_type, spatial_dimension, opt_param, *this, res_id); resolutions.push_back(std::move(resolution)); return *(resolutions.back()); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initResolutions() { AKANTU_DEBUG_ASSERT(resolutions.size() != 0, "No resolutions to initialize !"); if (!are_resolutions_instantiated) instantiateResolutions(); // \TODO check if each resolution needs a initResolution() method } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initModel() { - getFEEngine().initShapeFunctions(_not_ghost); - getFEEngine().initShapeFunctions(_ghost); + + AKANTU_DEBUG_IN(); + + getFEEngine("ContactMechanicsModel").initShapeFunctions(_not_ghost); + getFEEngine("ContactMechanicsModel").initShapeFunctions(_ghost); + + getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_not_ghost); + getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_ghost); + + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initSolver( TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType) { // 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->normal_force, spatial_dimension, "normal_force"); this->allocNodalField(this->tangential_force, spatial_dimension, "tangential_force"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->nodal_area, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); - this->allocNodalField(this->stick_or_slip, 1, "stick_or_slip"); + this->allocNodalField(this->contact_state, 1, "contact_state"); this->allocNodalField(this->previous_master_elements, 1, "previous_master_elements"); this->allocNodalField(this->normals, spatial_dimension, "normals"); auto surface_dimension = spatial_dimension - 1; this->allocNodalField(this->tangents, surface_dimension*spatial_dimension, "tangents"); this->allocNodalField(this->previous_tangents, surface_dimension*spatial_dimension, "previous_tangents"); - this->allocNodalField(this->projections, surface_dimension, - "projections"); - this->allocNodalField(this->previous_projections, surface_dimension, - "previous_projections"); - this->allocNodalField(this->stick_projections, surface_dimension, - "stick_projections"); this->allocNodalField(this->tangential_tractions, surface_dimension, "tangential_tractions"); this->allocNodalField(this->previous_tangential_tractions, surface_dimension, "previous_tangential_tractions"); // todo register multipliers as dofs for lagrange multipliers } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_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"); UInt nb_nodes = mesh.getNbNodes(); this->internal_force->clear(); this->normal_force->clear(); this->tangential_force->clear(); internal_force->resize(nb_nodes, 0.); normal_force->resize(nb_nodes, 0.); tangential_force->resize(nb_nodes, 0.); // 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); // TODO : uncomment when developing code for parallelization, // currently it addes the force twice for not ghost elements // hence source of error AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { // save the previous state this->savePreviousState(); contact_elements.clear(); contact_elements.resize(0); // this needs to be resized if cohesive elements are added UInt nb_nodes = mesh.getNbNodes(); auto resize_arrays = [&](auto & internal_array) { internal_array->clear(); internal_array->resize(nb_nodes, 0.); }; resize_arrays(gaps); resize_arrays(normals); resize_arrays(tangents); resize_arrays(projections); resize_arrays(tangential_tractions); - resize_arrays(stick_or_slip); + resize_arrays(contact_state); this->detector->search(contact_elements, *gaps, *normals, *projections); for (auto & gap : *gaps) { if (gap < 0) gap = std::abs(gap); else gap = -gap; } if (contact_elements.size() != 0) { this->computeNodalAreas(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::savePreviousState() { AKANTU_DEBUG_IN(); // saving previous natural projections previous_projections->clear(); previous_projections->resize(projections->size(), 0.); (*previous_projections).copy(*projections); // saving previous tangents previous_tangents->clear(); previous_tangents->resize(tangents->size(), 0.); (*previous_tangents).copy(*tangents); // saving previous tangential traction previous_tangential_tractions->clear(); previous_tangential_tractions->resize(tangential_tractions->size(), 0.); (*previous_tangential_tractions).copy(*tangential_tractions); previous_master_elements->clear(); previous_master_elements->resize(projections->size()); for (auto & element : contact_elements) { (*previous_master_elements)[element.slave] = element.master; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::computeNodalAreas() { UInt nb_nodes = mesh.getNbNodes(); nodal_area->clear(); external_force->clear(); nodal_area->resize(nb_nodes, 0.); external_force->resize(nb_nodes, 0.); - - auto & fem_boundary = this->getFEEngineBoundary(); + + auto & fem_boundary = this->getFEEngineBoundary("ContactMechanicsModel"); + fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); + + /*getFEEngine("ContactFacetsFEEngine").initShapeFunctions(getPositions(), _not_ghost); + getFEEngine("ContactFacetsFEEngine").initShapeFunctions(getPositions(), _ghost); - // TODO find a better method to compute nodal area + auto & fem_boundary = this->getFEEngineBoundary("ContactFacetsFEEngine"); + + fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); + fem_boundary.computeNormalsOnIntegrationPoints(_ghost);*/ + switch (spatial_dimension) { case 1: { for (auto && area : *nodal_area) { area = 1.; } break; } case 2: case 3: { this->applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), mesh.getElementGroup("contact_surface")); 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); } break; } default: break; } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ 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(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() { for (auto & resolution : resolutions) resolution->beforeSolveStep(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep(bool converged) { for (auto & resolution : resolutions) resolution->afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr 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["normal_force"] = this->normal_force; real_nodal_fields["tangential_force"] = this->tangential_force; real_nodal_fields["blocked_dofs"] = this->blocked_dofs; real_nodal_fields["normals"] = this->normals; real_nodal_fields["tangents"] = this->tangents; real_nodal_fields["gaps"] = this->gaps; real_nodal_fields["areas"] = this->nodal_area; - real_nodal_fields["stick_or_slip"] = this->stick_or_slip; + real_nodal_fields["contact_state"] = this->contact_state; std::shared_ptr field; 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 /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldReal(const std::string & , const std::string & , bool) { 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/contact_mechanics_model.hh b/src/model/contact_mechanics/contact_mechanics_model.hh index 08e1c6c5c..26677e0d1 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,386 +1,383 @@ /** * @file contact_mechanics_model.hh * * @author Mohit Pundir * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Contact Mechanics * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "contact_detector.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_MECHANICS_MODEL_HH__ #define __AKANTU_CONTACT_MECHANICS_MODEL_HH__ namespace akantu { class Resolution; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { + + + /* -------------------------------------------------------------------------- */ class ContactMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: ContactMechanicsModel( Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0, std::shared_ptr dof_manager = nullptr, const ModelType model_type = ModelType::_contact_mechanics_model); ~ContactMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// allocate all vectors void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize all internal arrays for resolutions void initResolutions(); /// initialize the modelType void initModel() override; /// call back for the solver, computes the force residual void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converted = true) override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Contact Detection */ /* ------------------------------------------------------------------------ */ public: void search(); void computeNodalAreas(); /* ------------------------------------------------------------------------ */ /* Contact Resolution */ /* ------------------------------------------------------------------------ */ public: /// register an empty contact resolution of a given type Resolution & registerNewResolution(const ID & res_name, const ID & res_type, const ID & opt_param); protected: /// register a resolution in the dynamic database Resolution & registerNewResolution(const ParserSection & res_section); /// read the resolution files to instantiate all the resolutions void instantiateResolutions(); /// save the parameters from previous state void savePreviousState(); /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: /// contact detection class friend class ContactDetector; /// contact resolution class friend class Resolution; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// get the ContactMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the ContactMechanicsModel::increment vector \warn only consistent /// if ContactMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the ContactMechanics::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the ContactMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the ContactMechanics::normal_force vector (normal forces) AKANTU_GET_MACRO(NormalForce, *normal_force, Array &); /// get the ContactMechanics::tangential_force vector (friction forces) AKANTU_GET_MACRO(TangentialForce, *tangential_force, Array &); /// get the ContactMechanics::traction vector (friction traction) AKANTU_GET_MACRO(TangentialTractions, *tangential_tractions, Array &); /// get the ContactMechanics::previous_tangential_tractions vector AKANTU_GET_MACRO(PreviousTangentialTractions, *previous_tangential_tractions, Array &); /// get the ContactMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the ContactMechanics::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the ContactMechanics::gaps (contact gaps) AKANTU_GET_MACRO(Gaps, *gaps, Array &); /// get the ContactMechanics::normals (normals on slave nodes) AKANTU_GET_MACRO(Normals, *normals, Array &); /// get the ContactMechanics::tangents (tangents on slave nodes) AKANTU_GET_MACRO(Tangents, *tangents, Array &); /// get the ContactMechanics::previous_tangents (tangents on slave nodes) AKANTU_GET_MACRO(PreviousTangents, *previous_tangents, Array &); /// get the ContactMechanics::areas (nodal areas) AKANTU_GET_MACRO(NodalArea, *nodal_area, Array &); - /// get the ContactMechanics::stick_projections (stick_projections) - AKANTU_GET_MACRO(StickProjections, *stick_projections, Array &); - /// get the ContactMechanics::previous_projections (previous_projections) AKANTU_GET_MACRO(PreviousProjections, *previous_projections, Array &); /// get the ContactMechanics::projections (projections) AKANTU_GET_MACRO(Projections, *projections, Array &); - /// get the ContactMechanics::stick_or_slip vector (slip/stick + /// get the ContactMechanics::contact_state vector (no_contact/stick/slip /// state) - AKANTU_GET_MACRO(StickSlip, *stick_or_slip, Array &); + AKANTU_GET_MACRO(ContactState, *contact_state, Array &); /// get the ContactMechanics::previous_master_elements AKANTU_GET_MACRO(PreviousMasterElements, *previous_master_elements, Array &); /// get contact detector AKANTU_GET_MACRO_NOT_CONST(ContactDetector, *detector, ContactDetector &); /// get the contact elements inline Array & getContactElements() { return contact_elements; } /// get the current positions of the nodes inline Array & getPositions() { return detector->getPositions(); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// tells if the resolutions are instantiated bool are_resolutions_instantiated; /// displacements array Array * displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// contact forces array Array * internal_force{nullptr}; /// external forces array Array * external_force{nullptr}; /// normal force array Array * normal_force{nullptr}; /// friction force array Array * tangential_force{nullptr}; /// friction traction array Array * tangential_tractions{nullptr}; /// previous friction traction array Array * previous_tangential_tractions{nullptr}; /// boundary vector Array * blocked_dofs{nullptr}; /// array to store gap between slave and master Array * gaps{nullptr}; /// array to store normals from master to slave Array * normals{nullptr}; /// array to store tangents on the master element Array * tangents{nullptr}; /// array to store previous tangents on the master element Array * previous_tangents{nullptr}; /// array to store nodal areas Array * nodal_area{nullptr}; /// array to store stick/slip state : - Array * stick_or_slip{nullptr}; - - /// array to store stick point projection in covariant basis - Array * stick_projections{nullptr}; + Array * contact_state{nullptr}; /// array to store previous projections in covariant basis Array * previous_projections{nullptr}; // array to store projections in covariant basis Array * projections{nullptr}; /// contact detection std::unique_ptr detector; /// list of contact resolutions std::vector> resolutions; /// mapping between resolution name and resolution internal id std::map resolutions_names_to_id; /// array to store contact elements Array contact_elements; /// array to store previous master elements Array * previous_master_elements{nullptr}; }; } // namespace akantu /* ------------------------------------------------------------------------ */ /* inline functions */ /* ------------------------------------------------------------------------ */ #include "parser.hh" #include "resolution.hh" /* ------------------------------------------------------------------------ */ #endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc index 060526984..50bad5424 100644 --- a/src/model/contact_mechanics/resolution.cc +++ b/src/model/contact_mechanics/resolution.cc @@ -1,252 +1,216 @@ /** * @file resolution.cc * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Implementation of common part of the contact resolution class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "resolution.hh" #include "contact_mechanics_model.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Resolution::Resolution(ContactMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_contact_resolution, id), fem(model.getFEEngine()), name(""), model(model) { AKANTU_DEBUG_IN(); spatial_dimension = model.getMesh().getSpatialDimension(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Resolution::~Resolution() = default; /* -------------------------------------------------------------------------- */ void Resolution::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("mu", mu, Real(0.), _pat_parsable | _pat_modifiable, "Friction Coefficient"); registerParam("is_master_deformable", is_master_deformable, bool(false), _pat_parsable | _pat_readable, "Is master surface deformable"); } /* -------------------------------------------------------------------------- */ void Resolution::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "Contact Resolution " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Resolution::assembleInternalForces(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); this->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleInternalForces() { AKANTU_DEBUG_IN(); for (auto & element : model.getContactElements()) { auto nb_nodes = element.getNbNodes(); Vector local_fn(nb_nodes * spatial_dimension); computeNormalForce(element, local_fn); Vector local_ft(nb_nodes * spatial_dimension); computeTangentialForce(element, local_ft); Vector local_fc(nb_nodes * spatial_dimension); local_fc = local_fn + local_ft; assembleLocalToGlobalArray(element, local_fn, model.getNormalForce()); assembleLocalToGlobalArray(element, local_ft, model.getTangentialForce()); assembleLocalToGlobalArray(element, local_fc, model.getInternalForce()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleLocalToGlobalArray(const ContactElement & element, Vector & local, Array & global) { auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn = const_cast(model.getMesh()).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; }; auto connectivity = get_connectivity(element.slave, element.master); UInt nb_dofs = global.getNbComponent(); UInt nb_nodes = is_master_deformable ? connectivity.size() : 1; for (UInt i : arange(nb_nodes)) { UInt n = connectivity[i]; for (UInt j : arange(nb_dofs)) { UInt offset_node = n * nb_dofs + j; global[offset_node] += local[i * nb_dofs + j]; } } } /* -------------------------------------------------------------------------- */ void Resolution::assembleStiffnessMatrix(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); auto & global_stiffness = const_cast(model.getDOFManager().getMatrix("K")); - - auto & gaps = model.getGaps(); - auto & projections = model.getProjections(); - auto & normals = model.getNormals(); - auto & tangents = model.getTangents(); - - UInt surface_dimension = spatial_dimension - 1; for (auto & element : model.getContactElements()) { auto nb_nodes = element.getNbNodes(); Matrix local_kn(nb_nodes * spatial_dimension, nb_nodes * spatial_dimension); computeNormalModuli(element, local_kn); assembleLocalToGlobalMatrix(element, local_kn, global_stiffness); - - - /*Real gap(gaps.begin()[element.slave]); - Vector normal(normals.begin(spatial_dimension)[element.slave]); - Vector projection(projections.begin(surface_dimension)[element.slave]); - Matrix covariant_basis(tangents.begin(surface_dimension, spatial_dimension)[element.slave]); - - //Matrix covariant_basis(surface_dimension, spatial_dimension); - //GeometryUtils::covariantBasis(model.getMesh(), model.getContactDetector().getPositions(), - // element.master, normal, projection, covariant_basis); - - Vector delta_g(nb_nodes * spatial_dimension); - ResolutionUtils::firstVariationNormalGap(element, projection, normal, delta_g); - - Matrix curvature(spatial_dimension, - surface_dimension * surface_dimension); - GeometryUtils::curvature(model.getMesh(), model.getContactDetector().getPositions(), - element.master, projection, curvature); - - Matrix ddelta_g(nb_nodes * spatial_dimension, nb_nodes * spatial_dimension); - ResolutionUtils::secondVariationNormalGap(element, covariant_basis, curvature, - projection, normal, gap, ddelta_g); - - Matrix local_kn(nb_nodes * spatial_dimension, nb_nodes * spatial_dimension); - computeNormalModuli(element, ddelta_g, delta_g, local_kn); - assembleLocalToGlobalMatrix(element, local_kn, global_stiffness); - + Matrix local_kt(nb_nodes * spatial_dimension, nb_nodes * spatial_dimension); - if (mu != 0) - computeTangentialModuli(element, ddelta_g, delta_g, local_kt); - assembleLocalToGlobalMatrix(element, local_kt, global_stiffness);*/ + computeTangentialModuli(element, local_kt); + assembleLocalToGlobalMatrix(element, local_kt, global_stiffness); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleLocalToGlobalMatrix(const ContactElement & element, const Matrix & local, SparseMatrix & global) { auto get_connectivity = [&](auto & slave, auto & master) { Vector master_conn = const_cast(model.getMesh()).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; }; auto connectivity = get_connectivity(element.slave, element.master); auto nb_dofs = spatial_dimension; UInt nb_nodes = is_master_deformable ? connectivity.size() : 1; UInt total_nb_dofs = nb_dofs * nb_nodes; std::vector equations; for (UInt i : arange(connectivity.size())) { UInt conn = connectivity[i]; for (UInt j : arange(nb_dofs)) { equations.push_back(conn * nb_dofs + j); } } for (UInt i : arange(total_nb_dofs)) { UInt row = equations[i]; for (UInt j : arange(total_nb_dofs)) { UInt col = equations[j]; global.add(row, col, local(i, j)); } } } /* -------------------------------------------------------------------------- */ -void Resolution::beforeSolveStep() { -} +void Resolution::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ -void Resolution::afterSolveStep(bool converged) { -} +void Resolution::afterSolveStep(bool converged) {} } // namespace akantu diff --git a/src/model/contact_mechanics/resolution.hh b/src/model/contact_mechanics/resolution.hh index 275d322af..afb7032f8 100644 --- a/src/model/contact_mechanics/resolution.hh +++ b/src/model/contact_mechanics/resolution.hh @@ -1,249 +1,251 @@ /** * @file resolution.hh * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Mother class for all contact resolutions * * @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_factory.hh" #include "aka_memory.hh" #include "parsable.hh" #include "parser.hh" #include "fe_engine.hh" #include "contact_element.hh" #include "resolution_utils.hh" #include "geometry_utils.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RESOLUTION_HH__ #define __AKANTU_RESOLUTION_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class ContactMechanicsModel; } // namespace akantu namespace akantu { /** * Interface of all contact resolutions * Prerequisites for a new resolution * - inherit from this class * - implement the following methods: * \code * * virtual void computeNormalForce(); - * virtual void computeFrictionForce(); + * virtual void computeTangentialForce(); + * virtual void computeNormalModuli(); + * virtual void computeTangentialModuli(); * * \endcode * */ class Resolution : public Memory, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructor/Destructor */ /* ------------------------------------------------------------------------ */ public: /// instantiate contact resolution with defaults Resolution(ContactMechanicsModel & model, const ID & id = ""); /// Destructor ~Resolution() override; protected: void initialize(); /// computes coordinates of a given element void computeCoordinates(const Element & , Matrix &); /* ------------------------------------------------------------------------ */ /* Functions that resolutions should reimplement for force */ /* ------------------------------------------------------------------------ */ public: /// computes the force vector due to normal traction virtual void computeNormalForce(__attribute__((unused)) const ContactElement &, __attribute__((unused)) Vector &) { AKANTU_TO_IMPLEMENT(); } /// computes the tangential force vector due to frictional traction virtual void computeTangentialForce(__attribute__((unused)) const ContactElement &, __attribute__((unused)) Vector &) { AKANTU_TO_IMPLEMENT(); } protected: /// local computation of trial tangential traction due to friction virtual void computeTrialTangentialTraction(__attribute__((unused)) const ContactElement &, __attribute__((unused)) const Matrix &, __attribute__((unused)) Vector &) { AKANTU_TO_IMPLEMENT(); } /// local computation of tangential traction due to stick virtual void computeStickTangentialTraction(__attribute__((unused)) const ContactElement &, __attribute__((unused)) Vector &, __attribute__((unused)) Vector &) { AKANTU_TO_IMPLEMENT(); } /// local computation of tangential traction due to slip virtual void computeSlipTangentialTraction(__attribute__((unused)) const ContactElement &, __attribute__((unused)) const Matrix &, __attribute__((unused)) Vector &, __attribute__((unused)) Vector &) { AKANTU_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* Functions that resolutions should reimplement for stiffness */ /* ------------------------------------------------------------------------ */ public: /// compute the normal moduli due to normal traction virtual void computeNormalModuli(__attribute__((unused)) const ContactElement &, __attribute__((unused)) Matrix & ) { AKANTU_TO_IMPLEMENT(); } /// compute the tangent moduli due to tangential traction virtual void computeTangentialModuli(__attribute__((unused)) const ContactElement &, __attribute__((unused)) Matrix & ) { AKANTU_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// assemble the residual for this resolution void assembleInternalForces(GhostType ghost_type); /// assemble the stiffness matrix for this resolution void assembleStiffnessMatrix(GhostType ghost_type); private: /// assemble the residual for this resolution void assembleInternalForces(); /// assemble the local array to global array for a contact element void assembleLocalToGlobalArray(const ContactElement & , Vector & , Array & ); /// assemble the local stiffness to global stiffness for a contact element void assembleLocalToGlobalMatrix(const ContactElement &, const Matrix &, SparseMatrix &); public: virtual void beforeSolveStep(); virtual void afterSolveStep(bool converged = true); public: /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// friction coefficient : mu Real mu; /// spatial dimension UInt spatial_dimension; /// is master surface deformable bool is_master_deformable; - /// Link to the fem object in the model + /// Link to the fe engine object in the model FEEngine & fem; /// resolution name std::string name; /// model to which the resolution belong ContactMechanicsModel & model; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Resolution & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { using ResolutionFactory = Factory; /// macaulay bracket to convert positive gap to zero template T macaulay(T var) {return var < 0 ? 0 : var; } template T heaviside(T var) {return var < 0 ? 0 : 1.0; } } // namespace akantu #define INSTANTIATE_RESOLUTION_ONLY(res_name) \ class res_name #define RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name) \ [](UInt dim, const ID &, ContactMechanicsModel & model, \ const ID & id) -> std::unique_ptr { \ switch (dim) { \ case 1: \ return std::make_unique(model, id); \ case 2: \ return std::make_unique(model, id); \ case 3: \ return std::make_unique(model, id); \ default: \ AKANTU_EXCEPTION("The dimension " \ << dim << "is not a valid dimension for the contact resolution " \ << #id); \ } \ } #define INSTANTIATE_RESOLUTION(id, res_name) \ INSTANTIATE_RESOLUTION_ONLY(res_name); \ static bool resolution_is_alocated_##id[[gnu::unused]] = \ ResolutionFactory::getInstance().registerAllocator( \ #id, RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name)) #endif /* __AKANTU_RESOLUTION_HH__ */ diff --git a/src/model/contact_mechanics/resolution_utils.cc b/src/model/contact_mechanics/resolution_utils.cc index 252270b6c..bf7ffd9dc 100644 --- a/src/model/contact_mechanics/resolution_utils.cc +++ b/src/model/contact_mechanics/resolution_utils.cc @@ -1,445 +1,445 @@ /** * @file resolution_utils.cc * * @author Mohit Pundir * * @date creation: Mon Mmay 20 2019 * @date last modification: Mon May 20 2019 * * @brief Implementation of various utilities neede for resolution class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "resolution_utils.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void ResolutionUtils::firstVariationNormalGap(const ContactElement & element, const Vector & projection, const Vector & normal, Vector & delta_g) { delta_g.clear(); const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; Vector shapes(Mesh::getNbNodesPerElement(type)); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL for (UInt i : arange(spatial_dimension)) { delta_g[i] = normal[i]; for (UInt j : arange(shapes.size())) { delta_g[(1 + j) * spatial_dimension + i] = -shapes[j] * normal[i]; } } } /* -------------------------------------------------------------------------- */ void ResolutionUtils::secondVariationNormalGap(const ContactElement & element, const Matrix & covariant_basis, const Matrix & curvature, const Vector & projection, const Vector & normal, Real & gap, Matrix & ddelta_g) { const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; UInt nb_nodes = element.getNbNodes(); Array dnds_n(nb_nodes * spatial_dimension, surface_dimension); ResolutionUtils::computeNalpha(element, projection, normal, dnds_n); Array delta_xi(nb_nodes * spatial_dimension, surface_dimension); ResolutionUtils::firstVariationNaturalCoordinate(element, covariant_basis, projection, normal, gap, delta_xi); Matrix a_alpha_beta(surface_dimension, surface_dimension); ResolutionUtils::computeMetricTensor(a_alpha_beta, covariant_basis); a_alpha_beta = a_alpha_beta.inverse(); Matrix h_alpha_beta(surface_dimension, surface_dimension); ResolutionUtils::computeSecondMetricTensor(element, curvature, normal, h_alpha_beta); for (auto && values : zip(arange(surface_dimension), make_view(dnds_n, dnds_n.size()), make_view(delta_xi, delta_xi.size()))) { auto & alpha = std::get<0>(values); auto & dnds_n_alpha = std::get<1>(values); auto & delta_xi_alpha = std::get<2>(values); // term 1 from Numerical methods in contact mechanics : Vlad // Yastrebov eq 2.48 Matrix mat_n(dnds_n_alpha.storage(), dnds_n_alpha.size(), 1); Matrix mat_xi(delta_xi_alpha.storage(), delta_xi_alpha.size(), 1); Matrix tmp1(dnds_n_alpha.size(), dnds_n_alpha.size()); tmp1.mul(mat_n, mat_xi, -1); Matrix tmp2(dnds_n_alpha.size(), dnds_n_alpha.size()); tmp2.mul(mat_xi, mat_n, -1); Matrix term1(dnds_n_alpha.size(), dnds_n_alpha.size()); term1 = tmp1 + tmp2; // computing term 2 & term 3 from Numerical methods in contact // mechanics : Vlad Yastrebov eq 2.48 Matrix term2(delta_xi_alpha.size(), delta_xi_alpha.size()); Matrix term3(dnds_n_alpha.size(), dnds_n_alpha.size()); for (auto && values2 : zip(arange(surface_dimension), make_view(dnds_n, dnds_n.size()), make_view(delta_xi, delta_xi.size()))) { auto & beta = std::get<0>(values2); auto & dnds_n_beta = std::get<1>(values2); auto & delta_xi_beta = std::get<2>(values2); // term 2 Matrix mat_xi_beta(delta_xi_beta.storage(), delta_xi.size(), 1); Matrix tmp3(delta_xi_beta.size(), delta_xi_beta.size()); Real pre_factor = h_alpha_beta(alpha, beta); for (auto k : arange(surface_dimension)) { for (auto m : arange(surface_dimension)) { pre_factor -= gap * h_alpha_beta(alpha, k) * a_alpha_beta(k, m) * h_alpha_beta(m, beta); } } pre_factor *= -1.; tmp3.mul(mat_xi, mat_xi_beta, pre_factor); // term 3 Matrix mat_n_beta(dnds_n_beta.storage(), dnds_n_beta.size(), 1); Real factor = gap * a_alpha_beta(alpha, beta); Matrix tmp4(dnds_n_alpha.size(), dnds_n_alpha.size()); tmp4.mul(mat_n, mat_n_beta, factor); term3 += tmp4; } ddelta_g += term1 + term2 + term3; } } /* -------------------------------------------------------------------------- */ void ResolutionUtils::computeTalpha(const ContactElement & element, const Matrix & covariant_basis, const Vector & projection, Array & t_alpha) { t_alpha.clear(); const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; Vector shapes(Mesh::getNbNodesPerElement(type)); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL for (auto && values : zip(covariant_basis.transpose(), make_view(t_alpha, t_alpha.size()))) { auto & tangent_beta = std::get<0>(values); auto & t_beta = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { t_beta[i] = tangent_beta(i); for (UInt j : arange(shapes.size())) { t_beta[(1 + j) * spatial_dimension + i] = -shapes[j] * tangent_beta(i); } } } } /* -------------------------------------------------------------------------- */ void ResolutionUtils::computeNalpha(const ContactElement & element, const Vector & projection, const Vector & normal, Array & n_alpha) { n_alpha.clear(); const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; Matrix shape_derivatives(surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shape_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL for (auto && values : zip(shape_derivatives.transpose(), make_view(n_alpha, n_alpha.size()))) { auto & dnds = std::get<0>(values); auto & n_s = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { n_s[i] = 0; for (UInt j : arange(dnds.size())) { n_s[(1 + j) * spatial_dimension + i] = -dnds(j) * normal[i]; } } } } /* -------------------------------------------------------------------------- */ void ResolutionUtils::firstVariationNaturalCoordinate(const ContactElement & element, const Matrix & covariant_basis, const Vector & projection, - const Vector & normal, const Real & gap, + const Vector & normal, const Real & /*gap*/, Array & delta_xi) { delta_xi.clear(); const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis); auto nb_nodes = element.getNbNodes(); Array t_alpha(nb_nodes * spatial_dimension, surface_dimension); Array n_alpha(nb_nodes * spatial_dimension, surface_dimension); ResolutionUtils::computeTalpha(element, covariant_basis, projection, t_alpha); ResolutionUtils::computeNalpha(element, projection, normal, n_alpha); for (auto && entry : zip(arange(surface_dimension), make_view(delta_xi, delta_xi.size()))) { auto & alpha = std::get<0>(entry); auto & d_alpha = std::get<1>(entry); for (auto && values : zip(arange(surface_dimension), make_view(t_alpha, t_alpha.size()), make_view(n_alpha, n_alpha.size()))) { auto & beta = std::get<0>(values); auto & t_beta = std::get<1>(values); - auto & n_beta = std::get<2>(values); + //auto & n_beta = std::get<2>(values); //d_alpha += (t_beta + gap * n_beta) * m_alpha_beta(alpha, //beta); d_alpha += t_beta * inv_A(alpha, beta); } } } - -/* -------------------------------------------------------------------------- */ -void ResolutionUtils::computeTalphabeta(Array & t_alpha_beta, - ContactElement & element) { - /*t_alpha_beta.clear(); - - const auto & type = element.master.type; - auto surface_dimension = Mesh::getSpatialDimension(type); - - Matrix shape_derivatives(surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(element.projection, shape_derivatives) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); - #undef GET_SHAPE_DERIVATIVES_NATURAL*/ - - //auto t_alpha_size = t_alpha_beta.size() * surface_dimension; - - //auto & tangents = element.tangents; - /*for(auto && entry : - zip(tangents.transpose(), - make_view(t_alpha_beta, t_alpha_size))) { - - auto & tangent_s = std::get<0>(entry); - auto & t_alpha = std::get<1>(entry); - - for(auto && values : - zip(shape_derivatives.transpose(), - make_view(t_alpha, t_alpha_beta.size()))) { - - auto & dnds = std::get<0>(values); - auto & t_alpha_s = std::get<1>(values); - - for (UInt i : arange(spatial_dimension)) { - t_alpha_s[i] = 0; - for (UInt j : arange(dnds.size())) { - t_alpha_s[(1 + j) * spatial_dimension + i] = -dnds(j) * tangent_s(i); - } - } - } - }*/ -} - -/* -------------------------------------------------------------------------- */ -void ResolutionUtils::computeNalphabeta(Array & n_alpha_beta, - ContactElement & /*element*/) { - n_alpha_beta.clear(); - - /*const auto & type = element.master.type; - auto surface_dimension = Mesh::getSpatialDimension(type); - auto spatial_dimension = surface_dimension + 1; - - Matrix shape_second_derivatives(surface_dimension * surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDN2DS2(element.projection, shape_second_derivatives) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); -#undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL - - for(auto && entry : - zip(shape_second_derivatives.transpose(), - make_view(n_alpha_beta, n_alpha_beta.size()))) { - - auto & dn2ds2 = std::get<0>(entry); - auto & n_alpha_s = std::get<1>(entry); - - for (UInt i : arange(spatial_dimension)) { - n_alpha_s[i] = 0; - for (UInt j : arange(dn2ds2.size())) { - n_alpha_s[(1 + j) * spatial_dimension + i] = -dn2ds2(j) * element.normal[i]; - } - } - - }*/ -} - -/* -------------------------------------------------------------------------- */ -void ResolutionUtils::computePalpha(Array & p_alpha, - ContactElement & /*element*/) { - p_alpha.clear(); - - /*const auto & type = element.master.type; - auto surface_dimension = Mesh::getSpatialDimension(type); - auto spatial_dimension = surface_dimension + 1; - - Matrix shape_derivatives(surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(element.projection, shape_derivatives) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL - - auto normalized_traction = element.traction/element.traction.norm(); - - for(auto && entry : - zip(shape_derivatives.transpose(), - make_view(p_alpha, p_alpha.size()))) { - auto & dnds = std::get<0>(entry); - auto & p_s = std::get<1>(entry); - - for(UInt i : arange(spatial_dimension)) { - p_s[i] = 0; - for(UInt j : arange(dnds.size())){ - p_s[(1 + j) * spatial_dimension + i] = -dnds(j) * normalized_traction[i]; - } - } - }*/ -} - -/* -------------------------------------------------------------------------- */ -void ResolutionUtils::computeGalpha(Array & g_alpha, Array & t_alpha_beta, - Array & d_alpha, Matrix & phi, - ContactElement & element) { - - g_alpha.clear(); - - /*const auto & type = element.master.type; - auto surface_dimension = Mesh::getSpatialDimension(type); - - auto & tangents = element.tangents; - auto tangential_gap = element.projection - element.previous_projection; - - for(auto alpha : arange(surface_dimension)) { - auto & g_a = g_alpha[alpha]; - auto & tangent_alpha = tangents[alpha]; - - for(auto beta : arange(surface_dimension)) { - auto & t_a_b = t_alpha_beta(alpha, beta); - auto & t_b_a = t_alpha_beta(beta, alpha); - auto & gt_beta = tangential_gap[beta]; - auto & tangents_beta = tangents[beta]; - - for(auto gamma : arange(surface_dimension)) { - auto & d_gamma = d_alpha(gamma); - auto tmp = phi(beta, gamma) * tangent_alpha + phi(alpha,gamma) * tangents_beta; - - g_a += (-t_a_b - t_b_a + tmp * d_gamma)*gt_beta; - } - } - - }*/ - -} - +// +///* -------------------------------------------------------------------------- */ +//void ResolutionUtils::computeTalphabeta(Array & t_alpha_beta, +// ContactElement & element) { +// /*t_alpha_beta.clear(); +// +// const auto & type = element.master.type; +// auto surface_dimension = Mesh::getSpatialDimension(type); +// +// Matrix shape_derivatives(surface_dimension, +// Mesh::getNbNodesPerElement(type)); +// +//#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ +// ElementClass::computeDNDS(element.projection, shape_derivatives) +// AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); +// #undef GET_SHAPE_DERIVATIVES_NATURAL*/ +// +// //auto t_alpha_size = t_alpha_beta.size() * surface_dimension; +// +// //auto & tangents = element.tangents; +// /*for(auto && entry : +// zip(tangents.transpose(), +// make_view(t_alpha_beta, t_alpha_size))) { +// +// auto & tangent_s = std::get<0>(entry); +// auto & t_alpha = std::get<1>(entry); +// +// for(auto && values : +// zip(shape_derivatives.transpose(), +// make_view(t_alpha, t_alpha_beta.size()))) { +// +// auto & dnds = std::get<0>(values); +// auto & t_alpha_s = std::get<1>(values); +// +// for (UInt i : arange(spatial_dimension)) { +// t_alpha_s[i] = 0; +// for (UInt j : arange(dnds.size())) { +// t_alpha_s[(1 + j) * spatial_dimension + i] = -dnds(j) * tangent_s(i); +// } +// } +// } +// }*/ +//} +// +///* -------------------------------------------------------------------------- */ +//void ResolutionUtils::computeNalphabeta(Array & n_alpha_beta, +// ContactElement & /*element*/) { +// n_alpha_beta.clear(); +// +// /*const auto & type = element.master.type; +// auto surface_dimension = Mesh::getSpatialDimension(type); +// auto spatial_dimension = surface_dimension + 1; +// +// Matrix shape_second_derivatives(surface_dimension * surface_dimension, +// Mesh::getNbNodesPerElement(type)); +// +//#define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ +// ElementClass::computeDN2DS2(element.projection, shape_second_derivatives) +// AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); +//#undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL +// +// for(auto && entry : +// zip(shape_second_derivatives.transpose(), +// make_view(n_alpha_beta, n_alpha_beta.size()))) { +// +// auto & dn2ds2 = std::get<0>(entry); +// auto & n_alpha_s = std::get<1>(entry); +// +// for (UInt i : arange(spatial_dimension)) { +// n_alpha_s[i] = 0; +// for (UInt j : arange(dn2ds2.size())) { +// n_alpha_s[(1 + j) * spatial_dimension + i] = -dn2ds2(j) * element.normal[i]; +// } +// } +// +// }*/ +//} +// +///* -------------------------------------------------------------------------- */ +//void ResolutionUtils::computePalpha(Array & p_alpha, +// ContactElement & /*element*/) { +// p_alpha.clear(); +// +// /*const auto & type = element.master.type; +// auto surface_dimension = Mesh::getSpatialDimension(type); +// auto spatial_dimension = surface_dimension + 1; +// +// Matrix shape_derivatives(surface_dimension, +// Mesh::getNbNodesPerElement(type)); +// +//#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ +// ElementClass::computeDNDS(element.projection, shape_derivatives) +// AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); +//#undef GET_SHAPE_DERIVATIVES_NATURAL +// +// auto normalized_traction = element.traction/element.traction.norm(); +// +// for(auto && entry : +// zip(shape_derivatives.transpose(), +// make_view(p_alpha, p_alpha.size()))) { +// auto & dnds = std::get<0>(entry); +// auto & p_s = std::get<1>(entry); +// +// for(UInt i : arange(spatial_dimension)) { +// p_s[i] = 0; +// for(UInt j : arange(dnds.size())){ +// p_s[(1 + j) * spatial_dimension + i] = -dnds(j) * normalized_traction[i]; +// } +// } +// }*/ +//} +// +///* -------------------------------------------------------------------------- */ +//void ResolutionUtils::computeGalpha(Array & g_alpha, Array & t_alpha_beta, +// Array & d_alpha, Matrix & phi, +// ContactElement & element) { +// +// g_alpha.clear(); +// +// /*const auto & type = element.master.type; +// auto surface_dimension = Mesh::getSpatialDimension(type); +// +// auto & tangents = element.tangents; +// auto tangential_gap = element.projection - element.previous_projection; +// +// for(auto alpha : arange(surface_dimension)) { +// auto & g_a = g_alpha[alpha]; +// auto & tangent_alpha = tangents[alpha]; +// +// for(auto beta : arange(surface_dimension)) { +// auto & t_a_b = t_alpha_beta(alpha, beta); +// auto & t_b_a = t_alpha_beta(beta, alpha); +// auto & gt_beta = tangential_gap[beta]; +// auto & tangents_beta = tangents[beta]; +// +// for(auto gamma : arange(surface_dimension)) { +// auto & d_gamma = d_alpha(gamma); +// auto tmp = phi(beta, gamma) * tangent_alpha + phi(alpha,gamma) * tangents_beta; +// +// g_a += (-t_a_b - t_b_a + tmp * d_gamma)*gt_beta; +// } +// } +// +// }*/ +// +//} +// /* -------------------------------------------------------------------------- */ void ResolutionUtils::computeMetricTensor(Matrix & m_alpha_beta, const Matrix & tangents) { m_alpha_beta.mul(tangents, tangents); } /* -------------------------------------------------------------------------- */ void ResolutionUtils::computeSecondMetricTensor(const ContactElement & element, const Matrix & curvature, const Vector & normal, Matrix & metric) { const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto i = 0; for (auto alpha : arange(surface_dimension) ) { for (auto beta : arange(surface_dimension)) { Vector temp(curvature(i)); metric(alpha, beta) = normal.dot(temp); i++; } } } } //akantu diff --git a/src/model/contact_mechanics/resolution_utils.hh b/src/model/contact_mechanics/resolution_utils.hh index 9869415da..c51c3b771 100644 --- a/src/model/contact_mechanics/resolution_utils.hh +++ b/src/model/contact_mechanics/resolution_utils.hh @@ -1,122 +1,122 @@ /** * @file resolution_utils.hh * * @author Mohit Pundir * * @date creation: Mon May 20 2019 * @date last modification: Mon May 20 2019 * * @brief All resolution utils necessary for various tasks * * @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 "contact_mechanics_model.hh" #include "contact_element.hh" #include "fe_engine.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RESOLUTION_UTILS_HH__ #define __AKANTU_RESOLUTION_UTILS_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class ResolutionUtils { /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// computes the first covariant metric tensor (@f$A_{\alpha\beta}@f$) where @f$\alpha, /// \beta@f$ are surface directions static void computeMetricTensor(Matrix & metric_tensor, const Matrix & tangents); /// computes the second covariant metric tensor /// (@f$H_{\alpha\beta}@f$) static void computeSecondMetricTensor(const ContactElement &, const Matrix &, const Vector &, Matrix &); /// computes the first variation of normal gap static void firstVariationNormalGap(const ContactElement & , const Vector &, const Vector &, Vector &); /// computes the seond variation of normal gap static void secondVariationNormalGap(const ContactElement & , const Matrix & , const Matrix &, const Vector &, const Vector &, Real &, Matrix & ); /// computes (@f$N_{\alpha}@f$) where \alpha is surface dimension /// and it is shape derivatives times normal static void computeNalpha(const ContactElement & , const Vector &, const Vector &, Array & ); /// computes (@f$T_{\alpha}@f$) where @f$\alpha@f$ is surface /// dimension and it is shape functions times the tangents static void computeTalpha(const ContactElement &, const Matrix &, const Vector &, Array & ); /// computes (@f$\nabla \xi_{\alpha}@f$) where @f$\alpha@f$ is surface /// dimension static void firstVariationNaturalCoordinate(const ContactElement &, const Matrix &, const Vector &, const Vector &, const Real &, Array &); ///computes second variation of surface parameter static void secondvariationNaturalCoordinate(const Vector & projection, const Vector & previous_projection, const Element & element, const Element & previous_element); /// computes @f$T_{\alpha\beta} @f$ which is shape derivatives /// times the tangents - static void computeTalphabeta(Array & t_alpha_beta, - ContactElement & element); + //static void computeTalphabeta(Array & t_alpha_beta, + // ContactElement & element); /// computes @f$N_{\alpha\beta} @f$ which is shape 2nd derivatives times /// the normal - static void computeNalphabeta(Array & n_alpha_beta, - ContactElement & element); + //static void computeNalphabeta(Array & n_alpha_beta, + // ContactElement & element); /// computes @f$P_{\alpha} @f$ - static void computePalpha(Array & p_alpha, ContactElement & element); + //static void computePalpha(Array & p_alpha, ContactElement & element); /// computes @f$G_{\alpha}@f$ - static void computeGalpha(Array & g_alpha, Array & t_alpha_beta, - Array & d_alpha, Matrix & phi, - ContactElement &); + //static void computeGalpha(Array & g_alpha, Array & t_alpha_beta, + // Array & d_alpha, Matrix & phi, + // ContactElement &); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class members */ /* ------------------------------------------------------------------------ */ private: }; } // namespace akantu #endif /* __AKANTU_RESOLUTION_UTILS_HH__ */ diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc index e1226d1fb..973c1e3f6 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc @@ -1,804 +1,865 @@ /** * @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_n", epsilon_n, Real(0.), _pat_parsable | _pat_modifiable, "Normal penalty parameter"); this->registerParam("epsilon_t", epsilon_t, Real(0.), _pat_parsable | _pat_modifiable, "Tangential penalty parameter"); } +/* -------------------------------------------------------------------------- */ +Real ResolutionPenalty::computeNormalTraction(Real & gap) { + return epsilon_n * macaulay(gap); +} + /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeNormalForce(const ContactElement & element, Vector & force) { force.clear(); auto & gaps = model.getGaps(); auto & projections = model.getProjections(); auto & normals = model.getNormals(); auto surface_dimension = spatial_dimension - 1; Real gap(gaps.begin()[element.slave]); Vector normal(normals.begin(spatial_dimension)[element.slave]); Vector projection(projections.begin(surface_dimension)[element.slave]); // compute normal traction - Real p_n = macaulay(gap) * epsilon_n; + Real p_n = computeNormalTraction(gap); // compute first variation of gap auto nb_nodes = element.getNbNodes(); Vector delta_gap(nb_nodes * spatial_dimension); ResolutionUtils::firstVariationNormalGap(element, projection, normal, delta_gap); // compute normal force auto & nodal_area = const_cast &>(model.getNodalArea()); for (UInt i : arange(force.size())) force[i] += delta_gap[i] * p_n * nodal_area[element.slave]; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTangentialForce(const ContactElement & element, Vector & force) { if (mu == 0) return; UInt surface_dimension = spatial_dimension - 1; // compute tangents auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); auto & tangents = model.getTangents(); Matrix covariant_basis(tangents.begin(surface_dimension, spatial_dimension)[element.slave]); GeometryUtils::covariantBasis(model.getMesh(), model.getContactDetector().getPositions(), element.master, normal, projection, covariant_basis); // check for no-contact to contact condition auto & previous_master_elements = model.getPreviousMasterElements(); auto & previous_element = previous_master_elements[element.slave]; if (previous_element.type == _not_defined) return; // compute tangential traction using return map algorithm auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction(tangential_tractions.begin(surface_dimension)[element.slave]); this->computeTangentialTraction(element, covariant_basis, tangential_traction); // compute first variation of natural coordinate auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; auto nb_nodes = element.getNbNodes(); Array delta_xi(nb_nodes * spatial_dimension, surface_dimension); ResolutionUtils::firstVariationNaturalCoordinate(element, covariant_basis, projection, normal, gap, delta_xi); // compute tangential force auto & nodal_area = const_cast &>(model.getNodalArea()); for (auto && values : zip(tangential_traction, make_view(delta_xi, delta_xi.size()))) { auto & traction_alpha = std::get<0>(values); auto & delta_xi_alpha = std::get<1>(values); force += delta_xi_alpha * traction_alpha * nodal_area[element.slave]; } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTangentialTraction(const ContactElement & element, const Matrix & covariant_basis, Vector & traction_tangential) { UInt surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; // Return map algorithm is employed // compute trial traction Vector traction_trial(surface_dimension); this->computeTrialTangentialTraction(element, covariant_basis, traction_trial); // compute norm of trial traction Real traction_trial_norm = 0; auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_trial_norm += traction_trial[i] * traction_trial[j] * inv_A(i, j); } } traction_trial_norm = sqrt(traction_trial_norm); // check stick or slip condition - auto & stick_or_slip = model.getStickSlip(); - auto & cond = stick_or_slip.begin()[element.slave]; + auto & contact_state = model.getContactState(); + auto & state = contact_state.begin()[element.slave]; - Real p_n = macaulay(gap) * epsilon_n; + Real p_n = computeNormalTraction(gap); bool stick = (traction_trial_norm <= mu * p_n) ? true : false; if (stick) { - cond = 1; + state = ContactState::_stick; computeStickTangentialTraction(element, traction_trial, traction_tangential); } else { - cond = 2; + state = ContactState::_slip; computeSlipTangentialTraction(element, covariant_basis, traction_trial, traction_tangential); } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTrialTangentialTraction(const ContactElement & element, const Matrix & covariant_basis, Vector & traction) { UInt surface_dimension = spatial_dimension - 1; - /*auto & projections = model.getProjections(); - Vector contravariant_projection(projections.begin(surface_dimension)[element.slave]); - - auto & stick_projections = model.getStickProjections(); - Vector covariant_stick(stick_projections.begin(surface_dimension)[element.slave]); - - Matrix contravariant_basis(surface_dimension, spatial_dimension); - GeometryUtils::contravariantBasis(covariant_basis, contravariant_basis); - - Vector covariant_projection(surface_dimension); - for (auto && values : zip(covariant_projection, - contravariant_projection, - contravariant_basis.transpose())) { - auto & temp = std::get<0>(values); - Vector contravariant_alpha(std::get<2>(values)); - - temp = contravariant_alpha.dot(contravariant_alpha); - temp *= std::get<1>(values); - } - - auto covariant_slip = covariant_projection - covariant_stick; - traction.mul(contravariant_basis, covariant_slip, epsilon_t);*/ - auto & projections = model.getProjections(); Vector current_projection(projections.begin(surface_dimension)[element.slave]); auto & previous_projections = model.getPreviousProjections(); Vector previous_projection(previous_projections.begin(surface_dimension)[element.slave]); // method from Laursen et. al. /*auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis); auto increment_projection = current_projection - previous_projection; traction.mul(covariant_metric_tensor, increment_projection, epsilon_t); auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); Vector previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]); traction = previous_traction + traction;*/ // method from Schweizerhof auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis); auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); Vector previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]); auto & previous_tangents = model.getPreviousTangents(); Matrix previous_covariant_basis(previous_tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto previous_contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(previous_covariant_basis); auto current_tangent = covariant_basis.transpose(); auto previous_tangent = previous_covariant_basis.transpose(); for (auto alpha :arange(surface_dimension)) { Vector tangent_alpha(current_tangent(alpha)); for (auto gamma : arange(surface_dimension)) { for (auto beta : arange(surface_dimension)) { Vector tangent_beta(previous_tangent(beta)); auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha); traction[alpha] += previous_traction[gamma]*previous_contravariant_metric_tensor(gamma, beta)*t_alpha_t_beta; } } } auto & previous_master_elements = model.getPreviousMasterElements(); auto & previous_element = previous_master_elements[element.slave]; Vector previous_real_projection(spatial_dimension); GeometryUtils::realProjection(model.getMesh(), model.getContactDetector().getPositions(), previous_element, previous_projection, previous_real_projection); Vector current_real_projection(spatial_dimension); GeometryUtils::realProjection(model.getMesh(), model.getContactDetector().getPositions(), element.master, current_projection, current_real_projection); auto increment_real = current_real_projection - previous_real_projection; Vector increment_xi(surface_dimension); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // increment in natural coordinate for (auto beta : arange(surface_dimension)) { for (auto gamma: arange(surface_dimension)) { auto temp = increment_real.dot(current_tangent(gamma)); temp *= contravariant_metric_tensor(beta, gamma); increment_xi[beta] += temp; } } Vector temp(surface_dimension); temp.mul(covariant_metric_tensor, increment_xi, epsilon_t); traction -= temp; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeStickTangentialTraction(const ContactElement & /*element*/, Vector & traction_trial, Vector & traction_tangential) { traction_tangential = traction_trial; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeSlipTangentialTraction(const ContactElement & element, const Matrix & covariant_basis, Vector & traction_trial, Vector & traction_tangential) { UInt surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; // compute norm of trial traction Real traction_trial_norm = 0; auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_trial_norm += traction_trial[i] * traction_trial[j] * inv_A(i, j); } } traction_trial_norm = sqrt(traction_trial_norm); auto slip_direction = traction_trial; slip_direction /= traction_trial_norm; - Real p_n = epsilon_n * macaulay(gap); + Real p_n = computeNormalTraction(gap); traction_tangential = slip_direction; - traction_tangential *= mu*p_n; - - /* - auto & stick_projections = model.getStickProjections(); - Vector covariant_stick(stick_projections.begin(surface_dimension)[element.slave]); - - - auto slip = macaulay(traction_trial_norm - mu * p_n); - slip /= epsilon_t; - - Vector covariant_slip(surface_dimension); - for (auto && values : zip(covariant_slip, - covariant_basis.transpose())) { - - auto & slip_alpha = std::get<0>(values); - Vector covariant_alpha(std::get<1>(values)); - slip_alpha = slip_direction.dot(covariant_alpha); - } - - covariant_stick += slip * covariant_slip;*/ + traction_tangential *= mu * p_n; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeNormalModuli(const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); Real gap(gaps.begin()[element.slave]); auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); auto & mesh = model.getMesh(); // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension*nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i*spatial_dimension + j) = 1; continue; } A(j, i*spatial_dimension + j) = -shapes[i-1]; } } // construct the main part of normal matrix Matrix k_main(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); Matrix n_outer_n(spatial_dimension, spatial_dimension); Matrix mat_n(normal.storage(), normal.size(), 1.); n_outer_n.mul(mat_n, mat_n); Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(n_outer_n, A); k_main.mul(A, tmp); k_main *= epsilon_n * heaviside(gap) * nodal_area; // construct the rotational part of the normal matrix auto & tangents = model.getTangents(); Matrix covariant_basis(tangents.begin(surface_dimension, spatial_dimension)[element.slave]); GeometryUtils::covariantBasis(model.getMesh(), model.getContactDetector().getPositions(), element.master, normal, projection, covariant_basis); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // computing shape derivatives Matrix shape_derivatives(surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shape_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL // consists of 2 rotational parts Matrix k_rot1(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); Matrix k_rot2(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); Matrix Aj(spatial_dimension, spatial_dimension*nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i*spatial_dimension + j) = 0; continue; } Aj(j, i*spatial_dimension + j) = dnds(i-1); } } }; for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent = std::get<1>(values1); Matrix n_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t(tangent.storage(), tangent.size(), 1.); n_outer_t.mul(mat_n, mat_t); Matrix t_outer_n(spatial_dimension, spatial_dimension); t_outer_n.mul(mat_t, mat_n); for (auto && values2 : enumerate(shape_derivatives.transpose())) { auto & beta = std::get<0>(values2); auto & dnds = std::get<1>(values2); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(n_outer_t, A); tmp1.mul(Aj, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_rot1 += tmp1; tmp.mul(t_outer_n, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_rot2 += tmp1; } } k_rot1 *= -epsilon_n * heaviside(gap) * gap * nodal_area; k_rot2 *= -epsilon_n * heaviside(gap) * gap * nodal_area; stiffness += k_main + k_rot1 + k_rot2; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTangentialModuli(const ContactElement & element, - Matrix & stiffness){ + Matrix & stiffness){ - - + if (mu == 0) + return; + + stiffness.clear(); + + auto & contact_state = model.getContactState(); + UInt state = contact_state.begin()[element.slave]; + + switch (state) { + case ContactState::_stick: { + computeStickModuli(element, stiffness); + break; + } + case ContactState::_slip: { + computeSlipModuli(element, stiffness); + break; + } + default: + break; + } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeStickModuli(const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; auto & mesh = model.getMesh(); // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension*nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i*spatial_dimension + j) = 1; continue; } A(j, i*spatial_dimension + j) = -shapes[i-1]; } } // computing shape derivatives Matrix shape_derivatives(surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shape_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Matrix Aj(spatial_dimension, spatial_dimension*nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i*spatial_dimension + j) = 0; continue; } Aj(j, i*spatial_dimension + j) = dnds(i-1); } } }; // tangents should have been calculated in normal modulii auto & tangents = model.getTangents(); Matrix covariant_basis(tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // construct 1st part of the stick modulii Matrix k_main(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); Matrix t_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); for (auto && values2 : enumerate(covariant_basis.transpose())) { auto & beta = std::get<0>(values2); auto & tangent_beta = std::get<1>(values2); Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); t_outer_t.mul(mat_t_alpha, mat_t_beta); Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(t_outer_t, A); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_main += tmp1; } } k_main *= -epsilon_t; // construct 2nd part of the stick modulii auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction(tangential_tractions.begin(surface_dimension)[element.slave]); Matrix k_second(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); for (auto alpha : arange(surface_dimension)) { Matrix k_sum(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); for (auto && values1 : enumerate(shape_derivatives.transpose())) { auto & beta = std::get<0>(values1); auto & dnds = std::get<1>(values1); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); for (auto && values2 : enumerate(covariant_basis.transpose())) { auto & gamma = std::get<0>(values2); auto & tangent_gamma = std::get<1>(values2); Matrix t_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), 1.); for (auto && values3 : enumerate(covariant_basis.transpose())) { auto & theta = std::get<0>(values3); auto & tangent_theta = std::get<1>(values3); Matrix mat_t_theta(tangent_theta.storage(), tangent_theta.size(), 1.); t_outer_t.mul(mat_t_gamma, mat_t_theta); Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(t_outer_t, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, theta) * contravariant_metric_tensor(beta, gamma); Matrix tmp2(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp3(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp2.mul(t_outer_t, A); tmp3.mul(Aj, tmp2); tmp3 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); k_sum += tmp1 + tmp3; } } } k_second += tangential_traction[alpha] * k_sum; } stiffness += k_main*nodal_area - k_second*nodal_area; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeSlipModuli(const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; + auto & gaps = model.getGaps(); + Real gap(gaps.begin()[element.slave]); + + auto & nodal_areas = model.getNodalArea(); + auto & nodal_area = nodal_areas.begin()[element.slave]; + + // compute normal traction + Real p_n = computeNormalTraction(gap); + auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); + // restructure normal as a matrix for an outer product Matrix mat_n(normal.storage(), normal.size(), 1.); auto & mesh = model.getMesh(); // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension*nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i*spatial_dimension + j) = 1; continue; } A(j, i*spatial_dimension + j) = -shapes[i-1]; } } // computing shape derivatives Matrix shape_derivatives(surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shape_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Matrix Aj(spatial_dimension, spatial_dimension*nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i*spatial_dimension + j) = 0; continue; } Aj(j, i*spatial_dimension + j) = dnds(i-1); } } }; // tangents should have been calculated in normal modulii auto & tangents = model.getTangents(); Matrix covariant_basis(tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction(tangential_tractions.begin(surface_dimension)[element.slave]); // compute norm of trial traction Real traction_norm = 0; auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_norm += tangential_traction[i] * tangential_traction[j] * inv_A(i, j); } } traction_norm = sqrt(traction_norm); - // construct 1st and 2nd part of stick modulii + // construct four parts of stick modulii (eq 107,107a-c) + Matrix k_first(nb_nodes_per_contact*spatial_dimension, + nb_nodes_per_contact*spatial_dimension); + Matrix k_second(nb_nodes_per_contact*spatial_dimension, + nb_nodes_per_contact*spatial_dimension); + Matrix k_third(nb_nodes_per_contact*spatial_dimension, + nb_nodes_per_contact*spatial_dimension); + Matrix k_fourth(nb_nodes_per_contact*spatial_dimension, + nb_nodes_per_contact*spatial_dimension); + + for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); + Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); + Matrix t_outer_n(spatial_dimension, spatial_dimension); - - for (auto && values2 : enumerate(covariant_basis.transpose())) { + Matrix t_outer_t(spatial_dimension, spatial_dimension); + + for (auto && values2 : zip(arange(surface_dimension), + covariant_basis.transpose(), + shape_derivatives.transpose())) { auto & beta = std::get<0>(values2); auto & tangent_beta = std::get<1>(values2); + auto & dnds = std::get<2>(values2); + //construct Aj from shape function wrt to jth natural + // coordinate + construct_Aj(dnds); + // eq 107 Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); t_outer_n.mul(mat_t_beta, mat_n); - + Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(t_outer_n, A); tmp1.mul(A, tmp); tmp1 *= epsilon_n * mu * tangential_traction[alpha] * contravariant_metric_tensor(alpha, beta); tmp1 /= traction_norm; - // add to t + k_first += tmp1 * nodal_area; + + // eq 107a + t_outer_t.mul(mat_t_alpha, mat_t_beta); + + tmp.mul(t_outer_t, A); + tmp1.mul(A, tmp); + + tmp1 *= epsilon_t * mu * p_n * contravariant_metric_tensor(alpha, beta); + tmp1 /= traction_norm; + + k_second += tmp1 * nodal_area; + for (auto && values3 : enumerate(covariant_basis.transpose())) { + auto & gamma = std::get<0>(values3); + auto & tangent_gamma = std::get<1>(values3); + + Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), 1.); + + for (auto && values4 : enumerate(covariant_basis.transpose())) { + auto & theta = std::get<0>(values4); + auto & tangent_theta = std::get<1>(values4); + + Matrix mat_t_theta(tangent_theta.storage(), tangent_theta.size(), 1.); + t_outer_t.mul(mat_t_gamma, mat_t_theta); + + // eq 107b + tmp.mul(t_outer_t, A); + tmp1.mul(A, tmp); + + tmp1 *= epsilon_t*mu*p_n*tangential_traction[alpha]*tangential_traction[beta]; + tmp1 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); + + k_third += tmp1 * nodal_area; + + // eq 107c + tmp.mul(t_outer_t, Aj); + tmp1.mul(A, tmp); + tmp1 *= contravariant_metric_tensor(alpha, theta) * contravariant_metric_tensor(beta, gamma); + tmp1 *= mu * p_n * tangential_traction[alpha]; + tmp1 /= traction_norm; + + Matrix tmp2(spatial_dimension, spatial_dimension*nb_nodes_per_contact); + Matrix tmp3(nb_nodes_per_contact*spatial_dimension, + spatial_dimension*nb_nodes_per_contact); + tmp2.mul(t_outer_t, A); + tmp3.mul(Aj, tmp2); + tmp3 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); + tmp3 *= mu * p_n * tangential_traction[alpha]; + tmp3 /= traction_norm; + + k_fourth += (tmp1 + tmp3) * nodal_area; + } + } } - - } - + stiffness += k_third + k_fourth - k_first - k_second; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::beforeSolveStep() { } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::afterSolveStep(bool converged) { /*auto method = model.getAnalysisMethod(); if (method == _explicit_lumped_mass) { return ; } auto & K = const_cast(model.getDOFManager().getMatrix("K")); auto k_min = K.min(); auto roundoff_error = 1e-17; const auto blocked_dofs = model.getDOFManager().getBlockedDOFs("displacement"); Real nb_unknowns = 0; for (auto & bld : make_view(blocked_dofs)) { if (not bld) nb_unknowns++; } auto max_epsilon_n = k_min / sqrt(nb_unknowns * roundoff_error); if (epsilon_n > max_epsilon_n) epsilon_n = max_epsilon_n;*/ } INSTANTIATE_RESOLUTION(penalty, ResolutionPenalty); } // namespace akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.hh b/src/model/contact_mechanics/resolutions/resolution_penalty.hh index 5e70fe17b..dfdbc00d4 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.hh +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.hh @@ -1,122 +1,125 @@ /** * @file contact_resolution_penalty.hh * * @author Mohit Pundir * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Jan 29 2018 * - * @brief Material isotropic thermo-elastic + * @brief Linear Penalty Resolution 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "resolution.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RESOLUTION_PENALTY_HH__ #define __AKANTU_RESOLUTION_PENALTY_HH__ namespace akantu { class ResolutionPenalty : public Resolution { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ResolutionPenalty(ContactMechanicsModel & model, const ID & id = ""); ~ResolutionPenalty() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize the resolution void initialize(); /// local computaion of stiffness matrix due to stick state void computeStickModuli(const ContactElement &, Matrix &); /// local computation of stiffness matrix due to slip state void computeSlipModuli(const ContactElement &, Matrix &); /* ------------------------------------------------------------------------ */ /* Methods for stiffness computation */ /* ------------------------------------------------------------------------ */ public: /// local computation of tangent moduli due to normal traction void computeNormalModuli(const ContactElement &, Matrix &) override; /// local computation of tangent moduli due to tangential traction void computeTangentialModuli(const ContactElement &, Matrix &) override; /* ------------------------------------------------------------------------ */ /* Methods for force computation */ /* ------------------------------------------------------------------------ */ public: /// local computation of normal force due to normal contact void computeNormalForce(const ContactElement &, Vector &) override; /// local computation of tangential force due to frictional traction void computeTangentialForce(const ContactElement &, Vector &) override; protected: + /// local computation of normal traction due to penetration + Real computeNormalTraction(Real &); + /// local computation of trial tangential traction due to friction void computeTrialTangentialTraction(const ContactElement &, const Matrix &, Vector &) override; /// local computation of tangential traction due to stick void computeStickTangentialTraction(const ContactElement &, Vector &, Vector &) override; /// local computation of tangential traction due to slip void computeSlipTangentialTraction(const ContactElement &, const Matrix &, Vector &, Vector &) override; /// local computation of tangential traction due to friction void computeTangentialTraction(const ContactElement &, const Matrix &, Vector &); public: void beforeSolveStep() override; void afterSolveStep(bool converged = true) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - /// penalty parameter for normal stress + /// penalty parameter for normal traction Real epsilon_n; - /// penalty parameter for tangential stress + /// penalty parameter for tangential traction Real epsilon_t; }; } // akantu #endif /* __AKANTU_RESOLUTION_PENALTY_HH__ */ diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc similarity index 81% copy from src/model/contact_mechanics/resolutions/resolution_penalty.cc copy to src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc index e1226d1fb..88cd48338 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc @@ -1,804 +1,863 @@ /** - * @file resolution_penalty.cc + * @file resolution_penalty_quadratic.cc * * @author Mohit Pundir * - * @date creation: Mon Jan 7 2019 - * @date last modification: Mon Jan 7 2019 + * @date creation: Sun Aug 2 2020 + * @date last modification: Sun Aug 2 2020 * - * @brief Specialization of the resolution class for the penalty method + * @brief Specialization of the resolution class for the quadratic 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" +#include "resolution_penalty_quadratic.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -ResolutionPenalty::ResolutionPenalty(ContactMechanicsModel & model, - const ID & id) - : Resolution(model, id) { +ResolutionPenaltyQuadratic::ResolutionPenaltyQuadratic(ContactMechanicsModel & model , + const ID & id) + : Parent(model, id) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::initialize() { - this->registerParam("epsilon_n", epsilon_n, 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 ResolutionPenaltyQuadratic::initialize() { + Parent::initialize(); } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeNormalForce(const ContactElement & element, - Vector & force) { - +Real ResolutionPenaltyQuadratic::computeNormalTraction(Real & gap) { + return epsilon_n *( macaulay(gap) * macaulay(gap) + macaulay(gap)); +} + +/* -------------------------------------------------------------------------- */ +void ResolutionPenaltyQuadratic::computeNormalForce(const ContactElement & element, + Vector & force) { force.clear(); - + auto & gaps = model.getGaps(); auto & projections = model.getProjections(); auto & normals = model.getNormals(); auto surface_dimension = spatial_dimension - 1; Real gap(gaps.begin()[element.slave]); Vector normal(normals.begin(spatial_dimension)[element.slave]); Vector projection(projections.begin(surface_dimension)[element.slave]); // compute normal traction - Real p_n = macaulay(gap) * epsilon_n; + Real p_n = computeNormalTraction(gap); // compute first variation of gap auto nb_nodes = element.getNbNodes(); Vector delta_gap(nb_nodes * spatial_dimension); ResolutionUtils::firstVariationNormalGap(element, projection, normal, delta_gap); // compute normal force auto & nodal_area = const_cast &>(model.getNodalArea()); for (UInt i : arange(force.size())) force[i] += delta_gap[i] * p_n * nodal_area[element.slave]; } - + /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentialForce(const ContactElement & element, - Vector & force) { +void ResolutionPenaltyQuadratic::computeTangentialForce(const ContactElement & element, + Vector & force) { if (mu == 0) return; + + force.clear(); UInt surface_dimension = spatial_dimension - 1; // compute tangents auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); auto & tangents = model.getTangents(); Matrix covariant_basis(tangents.begin(surface_dimension, spatial_dimension)[element.slave]); GeometryUtils::covariantBasis(model.getMesh(), model.getContactDetector().getPositions(), element.master, normal, projection, covariant_basis); // check for no-contact to contact condition auto & previous_master_elements = model.getPreviousMasterElements(); auto & previous_element = previous_master_elements[element.slave]; if (previous_element.type == _not_defined) return; // compute tangential traction using return map algorithm auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction(tangential_tractions.begin(surface_dimension)[element.slave]); this->computeTangentialTraction(element, covariant_basis, tangential_traction); // compute first variation of natural coordinate auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; auto nb_nodes = element.getNbNodes(); Array delta_xi(nb_nodes * spatial_dimension, surface_dimension); ResolutionUtils::firstVariationNaturalCoordinate(element, covariant_basis, projection, normal, gap, delta_xi); // compute tangential force auto & nodal_area = const_cast &>(model.getNodalArea()); for (auto && values : zip(tangential_traction, make_view(delta_xi, delta_xi.size()))) { auto & traction_alpha = std::get<0>(values); auto & delta_xi_alpha = std::get<1>(values); force += delta_xi_alpha * traction_alpha * nodal_area[element.slave]; } } + /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentialTraction(const ContactElement & element, +void ResolutionPenaltyQuadratic::computeTangentialTraction(const ContactElement & element, const Matrix & covariant_basis, Vector & traction_tangential) { UInt surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; // Return map algorithm is employed // compute trial traction Vector traction_trial(surface_dimension); this->computeTrialTangentialTraction(element, covariant_basis, traction_trial); // compute norm of trial traction Real traction_trial_norm = 0; auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_trial_norm += traction_trial[i] * traction_trial[j] * inv_A(i, j); } } traction_trial_norm = sqrt(traction_trial_norm); // check stick or slip condition - auto & stick_or_slip = model.getStickSlip(); - auto & cond = stick_or_slip.begin()[element.slave]; + auto & contact_state = model.getContactState(); + auto & state = contact_state.begin()[element.slave]; - Real p_n = macaulay(gap) * epsilon_n; + Real p_n = computeNormalTraction(gap); bool stick = (traction_trial_norm <= mu * p_n) ? true : false; if (stick) { - cond = 1; + state = ContactState::_stick; computeStickTangentialTraction(element, traction_trial, traction_tangential); } else { - cond = 2; + state = ContactState::_slip; computeSlipTangentialTraction(element, covariant_basis, traction_trial, traction_tangential); } } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTrialTangentialTraction(const ContactElement & element, +void ResolutionPenaltyQuadratic::computeTrialTangentialTraction(const ContactElement & element, const Matrix & covariant_basis, Vector & traction) { UInt surface_dimension = spatial_dimension - 1; - /*auto & projections = model.getProjections(); - Vector contravariant_projection(projections.begin(surface_dimension)[element.slave]); - - auto & stick_projections = model.getStickProjections(); - Vector covariant_stick(stick_projections.begin(surface_dimension)[element.slave]); - - Matrix contravariant_basis(surface_dimension, spatial_dimension); - GeometryUtils::contravariantBasis(covariant_basis, contravariant_basis); - - Vector covariant_projection(surface_dimension); - for (auto && values : zip(covariant_projection, - contravariant_projection, - contravariant_basis.transpose())) { - auto & temp = std::get<0>(values); - Vector contravariant_alpha(std::get<2>(values)); - - temp = contravariant_alpha.dot(contravariant_alpha); - temp *= std::get<1>(values); - } - - auto covariant_slip = covariant_projection - covariant_stick; - traction.mul(contravariant_basis, covariant_slip, epsilon_t);*/ - auto & projections = model.getProjections(); Vector current_projection(projections.begin(surface_dimension)[element.slave]); auto & previous_projections = model.getPreviousProjections(); Vector previous_projection(previous_projections.begin(surface_dimension)[element.slave]); // method from Laursen et. al. /*auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis); auto increment_projection = current_projection - previous_projection; traction.mul(covariant_metric_tensor, increment_projection, epsilon_t); auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); Vector previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]); traction = previous_traction + traction;*/ // method from Schweizerhof auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis); auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); Vector previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]); auto & previous_tangents = model.getPreviousTangents(); Matrix previous_covariant_basis(previous_tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto previous_contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(previous_covariant_basis); auto current_tangent = covariant_basis.transpose(); auto previous_tangent = previous_covariant_basis.transpose(); for (auto alpha :arange(surface_dimension)) { Vector tangent_alpha(current_tangent(alpha)); for (auto gamma : arange(surface_dimension)) { for (auto beta : arange(surface_dimension)) { Vector tangent_beta(previous_tangent(beta)); auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha); traction[alpha] += previous_traction[gamma]*previous_contravariant_metric_tensor(gamma, beta)*t_alpha_t_beta; } } } auto & previous_master_elements = model.getPreviousMasterElements(); auto & previous_element = previous_master_elements[element.slave]; Vector previous_real_projection(spatial_dimension); GeometryUtils::realProjection(model.getMesh(), model.getContactDetector().getPositions(), previous_element, previous_projection, previous_real_projection); Vector current_real_projection(spatial_dimension); GeometryUtils::realProjection(model.getMesh(), model.getContactDetector().getPositions(), element.master, current_projection, current_real_projection); auto increment_real = current_real_projection - previous_real_projection; Vector increment_xi(surface_dimension); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // increment in natural coordinate for (auto beta : arange(surface_dimension)) { for (auto gamma: arange(surface_dimension)) { auto temp = increment_real.dot(current_tangent(gamma)); temp *= contravariant_metric_tensor(beta, gamma); increment_xi[beta] += temp; } } Vector temp(surface_dimension); temp.mul(covariant_metric_tensor, increment_xi, epsilon_t); traction -= temp; } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeStickTangentialTraction(const ContactElement & /*element*/, +void ResolutionPenaltyQuadratic::computeStickTangentialTraction(const ContactElement & /*element*/, Vector & traction_trial, Vector & traction_tangential) { traction_tangential = traction_trial; } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeSlipTangentialTraction(const ContactElement & element, +void ResolutionPenaltyQuadratic::computeSlipTangentialTraction(const ContactElement & element, const Matrix & covariant_basis, Vector & traction_trial, Vector & traction_tangential) { UInt surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; // compute norm of trial traction Real traction_trial_norm = 0; auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_trial_norm += traction_trial[i] * traction_trial[j] * inv_A(i, j); } } traction_trial_norm = sqrt(traction_trial_norm); auto slip_direction = traction_trial; slip_direction /= traction_trial_norm; - Real p_n = epsilon_n * macaulay(gap); + Real p_n = computeNormalTraction(gap); traction_tangential = slip_direction; - traction_tangential *= mu*p_n; - - /* - auto & stick_projections = model.getStickProjections(); - Vector covariant_stick(stick_projections.begin(surface_dimension)[element.slave]); - - - auto slip = macaulay(traction_trial_norm - mu * p_n); - slip /= epsilon_t; - - Vector covariant_slip(surface_dimension); - for (auto && values : zip(covariant_slip, - covariant_basis.transpose())) { - - auto & slip_alpha = std::get<0>(values); - Vector covariant_alpha(std::get<1>(values)); - slip_alpha = slip_direction.dot(covariant_alpha); - } - - covariant_stick += slip * covariant_slip;*/ + traction_tangential *= mu * p_n; } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeNormalModuli(const ContactElement & element, +void ResolutionPenaltyQuadratic::computeNormalModuli(const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); Real gap(gaps.begin()[element.slave]); auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); auto & mesh = model.getMesh(); // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension*nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i*spatial_dimension + j) = 1; continue; } A(j, i*spatial_dimension + j) = -shapes[i-1]; } } // construct the main part of normal matrix Matrix k_main(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); Matrix n_outer_n(spatial_dimension, spatial_dimension); Matrix mat_n(normal.storage(), normal.size(), 1.); n_outer_n.mul(mat_n, mat_n); Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(n_outer_n, A); k_main.mul(A, tmp); k_main *= epsilon_n * heaviside(gap) * nodal_area; // construct the rotational part of the normal matrix auto & tangents = model.getTangents(); Matrix covariant_basis(tangents.begin(surface_dimension, spatial_dimension)[element.slave]); GeometryUtils::covariantBasis(model.getMesh(), model.getContactDetector().getPositions(), element.master, normal, projection, covariant_basis); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // computing shape derivatives Matrix shape_derivatives(surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shape_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL // consists of 2 rotational parts Matrix k_rot1(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); Matrix k_rot2(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); Matrix Aj(spatial_dimension, spatial_dimension*nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i*spatial_dimension + j) = 0; continue; } Aj(j, i*spatial_dimension + j) = dnds(i-1); } } }; for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent = std::get<1>(values1); Matrix n_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t(tangent.storage(), tangent.size(), 1.); n_outer_t.mul(mat_n, mat_t); Matrix t_outer_n(spatial_dimension, spatial_dimension); t_outer_n.mul(mat_t, mat_n); for (auto && values2 : enumerate(shape_derivatives.transpose())) { auto & beta = std::get<0>(values2); auto & dnds = std::get<1>(values2); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(n_outer_t, A); tmp1.mul(Aj, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_rot1 += tmp1; tmp.mul(t_outer_n, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_rot2 += tmp1; } } k_rot1 *= -epsilon_n * heaviside(gap) * gap * nodal_area; k_rot2 *= -epsilon_n * heaviside(gap) * gap * nodal_area; stiffness += k_main + k_rot1 + k_rot2; } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentialModuli(const ContactElement & element, - Matrix & stiffness){ +void ResolutionPenaltyQuadratic::computeTangentialModuli(const ContactElement & element, + Matrix & stiffness){ - - + if (mu == 0) + return; + + stiffness.clear(); + + auto & contact_state = model.getContactState(); + UInt state = contact_state.begin()[element.slave]; + + switch (state) { + case ContactState::_stick: { + computeStickModuli(element, stiffness); + break; + } + case ContactState::_slip: { + computeSlipModuli(element, stiffness); + break; + } + default: + break; + } } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeStickModuli(const ContactElement & element, +void ResolutionPenaltyQuadratic::computeStickModuli(const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; auto & mesh = model.getMesh(); // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension*nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i*spatial_dimension + j) = 1; continue; } A(j, i*spatial_dimension + j) = -shapes[i-1]; } } // computing shape derivatives Matrix shape_derivatives(surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shape_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Matrix Aj(spatial_dimension, spatial_dimension*nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i*spatial_dimension + j) = 0; continue; } Aj(j, i*spatial_dimension + j) = dnds(i-1); } } }; // tangents should have been calculated in normal modulii auto & tangents = model.getTangents(); Matrix covariant_basis(tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // construct 1st part of the stick modulii Matrix k_main(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); Matrix t_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); for (auto && values2 : enumerate(covariant_basis.transpose())) { auto & beta = std::get<0>(values2); auto & tangent_beta = std::get<1>(values2); Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); t_outer_t.mul(mat_t_alpha, mat_t_beta); Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(t_outer_t, A); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_main += tmp1; } } k_main *= -epsilon_t; // construct 2nd part of the stick modulii auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction(tangential_tractions.begin(surface_dimension)[element.slave]); Matrix k_second(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); for (auto alpha : arange(surface_dimension)) { Matrix k_sum(nb_nodes_per_contact*spatial_dimension, nb_nodes_per_contact*spatial_dimension); for (auto && values1 : enumerate(shape_derivatives.transpose())) { auto & beta = std::get<0>(values1); auto & dnds = std::get<1>(values1); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); for (auto && values2 : enumerate(covariant_basis.transpose())) { auto & gamma = std::get<0>(values2); auto & tangent_gamma = std::get<1>(values2); Matrix t_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), 1.); for (auto && values3 : enumerate(covariant_basis.transpose())) { auto & theta = std::get<0>(values3); auto & tangent_theta = std::get<1>(values3); Matrix mat_t_theta(tangent_theta.storage(), tangent_theta.size(), 1.); t_outer_t.mul(mat_t_gamma, mat_t_theta); Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(t_outer_t, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, theta) * contravariant_metric_tensor(beta, gamma); Matrix tmp2(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp3(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp2.mul(t_outer_t, A); tmp3.mul(Aj, tmp2); tmp3 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); k_sum += tmp1 + tmp3; } } } k_second += tangential_traction[alpha] * k_sum; } stiffness += k_main*nodal_area - k_second*nodal_area; } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeSlipModuli(const ContactElement & element, - Matrix & stiffness) { +void ResolutionPenaltyQuadratic::computeSlipModuli(const ContactElement & element, + Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; + auto & gaps = model.getGaps(); + Real gap(gaps.begin()[element.slave]); + + auto & nodal_areas = model.getNodalArea(); + auto & nodal_area = nodal_areas.begin()[element.slave]; + + // compute normal traction + Real p_n = computeNormalTraction(gap); + auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); + // restructure normal as a matrix for an outer product Matrix mat_n(normal.storage(), normal.size(), 1.); auto & mesh = model.getMesh(); // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension*nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i*spatial_dimension + j) = 1; continue; } A(j, i*spatial_dimension + j) = -shapes[i-1]; } } // computing shape derivatives Matrix shape_derivatives(surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shape_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Matrix Aj(spatial_dimension, spatial_dimension*nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i*spatial_dimension + j) = 0; continue; } Aj(j, i*spatial_dimension + j) = dnds(i-1); } } }; // tangents should have been calculated in normal modulii auto & tangents = model.getTangents(); Matrix covariant_basis(tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction(tangential_tractions.begin(surface_dimension)[element.slave]); // compute norm of trial traction Real traction_norm = 0; auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_norm += tangential_traction[i] * tangential_traction[j] * inv_A(i, j); } } traction_norm = sqrt(traction_norm); - // construct 1st and 2nd part of stick modulii + // construct four parts of stick modulii (eq 107,107a-c) + Matrix k_first(nb_nodes_per_contact*spatial_dimension, + nb_nodes_per_contact*spatial_dimension); + Matrix k_second(nb_nodes_per_contact*spatial_dimension, + nb_nodes_per_contact*spatial_dimension); + Matrix k_third(nb_nodes_per_contact*spatial_dimension, + nb_nodes_per_contact*spatial_dimension); + Matrix k_fourth(nb_nodes_per_contact*spatial_dimension, + nb_nodes_per_contact*spatial_dimension); + + for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); + Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); + Matrix t_outer_n(spatial_dimension, spatial_dimension); - - for (auto && values2 : enumerate(covariant_basis.transpose())) { + Matrix t_outer_t(spatial_dimension, spatial_dimension); + + for (auto && values2 : zip(arange(surface_dimension), + covariant_basis.transpose(), + shape_derivatives.transpose())) { auto & beta = std::get<0>(values2); auto & tangent_beta = std::get<1>(values2); + auto & dnds = std::get<2>(values2); + //construct Aj from shape function wrt to jth natural + // coordinate + construct_Aj(dnds); + // eq 107 Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); t_outer_n.mul(mat_t_beta, mat_n); - + Matrix tmp(spatial_dimension, spatial_dimension*nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact*spatial_dimension, spatial_dimension*nb_nodes_per_contact); tmp.mul(t_outer_n, A); tmp1.mul(A, tmp); tmp1 *= epsilon_n * mu * tangential_traction[alpha] * contravariant_metric_tensor(alpha, beta); tmp1 /= traction_norm; - // add to t + k_first += tmp1 * nodal_area; + + // eq 107a + t_outer_t.mul(mat_t_alpha, mat_t_beta); + tmp.mul(t_outer_t, A); + tmp1.mul(A, tmp); + + tmp1 *= epsilon_t * mu * p_n * contravariant_metric_tensor(alpha, beta); + tmp1 /= traction_norm; + + k_second += tmp1 * nodal_area; + + for (auto && values3 : enumerate(covariant_basis.transpose())) { + auto & gamma = std::get<0>(values3); + auto & tangent_gamma = std::get<1>(values3); + + Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), 1.); + + for (auto && values4 : enumerate(covariant_basis.transpose())) { + auto & theta = std::get<0>(values4); + auto & tangent_theta = std::get<1>(values4); + + Matrix mat_t_theta(tangent_theta.storage(), tangent_theta.size(), 1.); + t_outer_t.mul(mat_t_gamma, mat_t_theta); + + // eq 107b + tmp.mul(t_outer_t, A); + tmp1.mul(A, tmp); + + tmp1 *= epsilon_t*mu*p_n*tangential_traction[alpha]*tangential_traction[beta]; + tmp1 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); + + k_third += tmp1 * nodal_area; + + // eq 107c + tmp.mul(t_outer_t, Aj); + tmp1.mul(A, tmp); + tmp1 *= contravariant_metric_tensor(alpha, theta) * contravariant_metric_tensor(beta, gamma); + tmp1 *= mu * p_n * tangential_traction[alpha]; + tmp1 /= traction_norm; + + Matrix tmp2(spatial_dimension, spatial_dimension*nb_nodes_per_contact); + Matrix tmp3(nb_nodes_per_contact*spatial_dimension, + spatial_dimension*nb_nodes_per_contact); + tmp2.mul(t_outer_t, A); + tmp3.mul(Aj, tmp2); + tmp3 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); + tmp3 *= mu * p_n * tangential_traction[alpha]; + tmp3 /= traction_norm; + + k_fourth += (tmp1 + tmp3) * nodal_area; + } + } } - - } - + stiffness += k_third + k_fourth - k_first - k_second; } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::beforeSolveStep() { +void ResolutionPenaltyQuadratic::beforeSolveStep() { } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::afterSolveStep(bool converged) { +void ResolutionPenaltyQuadratic::afterSolveStep(bool converged) { /*auto method = model.getAnalysisMethod(); if (method == _explicit_lumped_mass) { return ; } auto & K = const_cast(model.getDOFManager().getMatrix("K")); auto k_min = K.min(); auto roundoff_error = 1e-17; const auto blocked_dofs = model.getDOFManager().getBlockedDOFs("displacement"); Real nb_unknowns = 0; for (auto & bld : make_view(blocked_dofs)) { if (not bld) nb_unknowns++; } auto max_epsilon_n = k_min / sqrt(nb_unknowns * roundoff_error); if (epsilon_n > max_epsilon_n) epsilon_n = max_epsilon_n;*/ } -INSTANTIATE_RESOLUTION(penalty, ResolutionPenalty); +INSTANTIATE_RESOLUTION(penalty_quadratic, ResolutionPenaltyQuadratic); + } // namespace akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.hh b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh similarity index 80% copy from src/model/contact_mechanics/resolutions/resolution_penalty.hh copy to src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh index 5e70fe17b..e02b38433 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.hh +++ b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.hh @@ -1,122 +1,119 @@ /** - * @file contact_resolution_penalty.hh + * @file contact_resolution_penalty_quadratic.hh * * @author Mohit Pundir * - * @date creation: Fri Jun 18 2010 - * @date last modification: Mon Jan 29 2018 + * @date creation: Sun Aug 02 2020 + * @date last modification: Sun Aug 02 2020 * - * @brief Material isotropic thermo-elastic + * @brief Quadratic Penalty Resolution 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" -#include "resolution.hh" +#include "resolution_penalty.hh" /* -------------------------------------------------------------------------- */ -#ifndef __AKANTU_RESOLUTION_PENALTY_HH__ -#define __AKANTU_RESOLUTION_PENALTY_HH__ +#ifndef __AKANTU_RESOLUTION_PENALTY_QUADRATIC_HH__ +#define __AKANTU_RESOLUTION_PENALTY_QUADRATIC_HH__ namespace akantu { -class ResolutionPenalty : public Resolution { + +class ResolutionPenaltyQuadratic : public ResolutionPenalty { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ +private: + using Parent = ResolutionPenalty; + public: - ResolutionPenalty(ContactMechanicsModel & model, const ID & id = ""); + ResolutionPenaltyQuadratic(ContactMechanicsModel & model, const ID & id = ""); - ~ResolutionPenalty() override = default; + ~ResolutionPenaltyQuadratic() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize the resolution void initialize(); /// local computaion of stiffness matrix due to stick state void computeStickModuli(const ContactElement &, Matrix &); /// local computation of stiffness matrix due to slip state void computeSlipModuli(const ContactElement &, Matrix &); /* ------------------------------------------------------------------------ */ /* Methods for stiffness computation */ /* ------------------------------------------------------------------------ */ public: /// local computation of tangent moduli due to normal traction void computeNormalModuli(const ContactElement &, Matrix &) override; /// local computation of tangent moduli due to tangential traction void computeTangentialModuli(const ContactElement &, Matrix &) override; /* ------------------------------------------------------------------------ */ /* Methods for force computation */ /* ------------------------------------------------------------------------ */ public: /// local computation of normal force due to normal contact void computeNormalForce(const ContactElement &, Vector &) override; /// local computation of tangential force due to frictional traction void computeTangentialForce(const ContactElement &, Vector &) override; protected: + /// local computation of normal traction due to penetration + Real computeNormalTraction(Real &); + /// local computation of trial tangential traction due to friction void computeTrialTangentialTraction(const ContactElement &, const Matrix &, Vector &) override; /// local computation of tangential traction due to stick void computeStickTangentialTraction(const ContactElement &, Vector &, Vector &) override; /// local computation of tangential traction due to slip void computeSlipTangentialTraction(const ContactElement &, const Matrix &, Vector &, Vector &) override; /// local computation of tangential traction due to friction void computeTangentialTraction(const ContactElement &, const Matrix &, Vector &); public: void beforeSolveStep() override; - void afterSolveStep(bool converged = true) override; - - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ -protected: - /// penalty parameter for normal stress - Real epsilon_n; - - /// penalty parameter for tangential stress - Real epsilon_t; + void afterSolveStep(bool converged = true) override; }; - } // akantu -#endif /* __AKANTU_RESOLUTION_PENALTY_HH__ */ + +#endif /* __AKANTU_RESOLUTION_PENALTY_QUADRATIC_HH__ */