diff --git a/doc/dev-doc/manual/contactmechanicsinternodesmodel.rst b/doc/dev-doc/manual/contactmechanicsinternodesmodel.rst new file mode 100644 index 000000000..0393c29d7 --- /dev/null +++ b/doc/dev-doc/manual/contactmechanicsinternodesmodel.rst @@ -0,0 +1,68 @@ +Contact Mechanics Internodes Model +================================== + +The contact mechanics internodes model is a implementation of +:cpp:class:`Model ` interface to handle contact between +two solid bodies using the INTERNODES (reference) method. + +Overview INTERNODES method +-------------------- + +The INTERNODES method in the context of cantact mechanics is algorithm to solve the +contact problem using nodes belonging to the interface. Information such as +displacement is transfered from one body to the other with a interpolation +using rescaled radial basis function. + +The initial configuration should already exhibit some interprenation to start +the algorithm. A set of predefined possible contact nodes of both bodies will +be taken in consideration to solve the problem. The algorithm will solve the +intenodes method. + + +Using the contact mechanics internodes model +-------------------------------------------- + +The :cpp:class:`ContactMechanicsInternodesModel +` works simalarly to the +:cpp:class:`ContactMechanicssModel `. +The class can be instanciated with an existing (see \ref{sect:common:mesh}) as follows:: + + ContactMechanicsInternodesModel contact_model(mesh, spatial_dimension); + +In contrast to the the standard :cpp:class:`ContactMechanicssModel +` the INTERNODES model handles the +contact resolution with the :cpp:class:`ContactDetectorInternodes +` and the coupling with a +solid model. Currently this method is fixed to use +:cpp:class:`SolidMechanicsModel `. + + +The initilization can be done with:: + + contact_model.initFull(_analysis_method =_static); + +The code belows gives the associated detector and solid model:: + + detector = contact_model.getContactDetector() + solid = contact_model.getSolidMechanicsModel() + +The boundary conditions can then be set directly through the solid model. + + +To initizial the contact detection of the model a configuration section can be +added to the 'material.dat' file. The possible nodes for the two interfaces can +be specified as follows. + +.. code-block:: + + contact_detector [ + master = contact_bottom + slave = contact_top + ] + + +Currently, one iteration of the INTERNODES method is possible:: + +.. code-block:: c++ + + contact_model.solveStep(); diff --git a/test/test_model/test_contact_mechanics_internodes_model/CMakeLists.txt b/test/test_model/test_contact_mechanics_internodes_model/CMakeLists.txt new file mode 100644 index 000000000..f7e5186e7 --- /dev/null +++ b/test/test_model/test_contact_mechanics_internodes_model/CMakeLists.txt @@ -0,0 +1,40 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Moritz Waldleben +# +# @date creation: Thu 16 2022 +# @date last modification: Thue 16 2022 +# +# @brief configuration python contact mechanics internodes test +# +# +# @section LICENSE +# +# Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# Akantu is free software: you can redistribute it and/or modify it under the +# terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License along +# with Akantu. If not, see . +# +#=============================================================================== + +add_mesh(contact_mesh prototype_internodes/contact.geo DIM 2 ORDER 1 OUTPUT contact.msh) + +register_test(test_internodes + SCRIPT test_internodes.py + PYTHON + FILES_TO_COPY prototype_internodes/material.dat + DEPENDS contact_mesh + PACKAGE python_interface contact_mechanics_internodes_model + ) diff --git a/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/example_direct.py b/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/example_direct.py index b9c97675a..d808256aa 100644 --- a/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/example_direct.py +++ b/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/example_direct.py @@ -1,153 +1,154 @@ #!/usr/bin/env python # coding: utf-8 import akantu as aka import numpy as np import scipy as sp import matplotlib.pyplot as plt import sys sys.path.append("..") from prototype_internodes.functions import * from prototype_internodes.functions_contact_probl import * from prototype_internodes.init_model import init_model np.set_printoptions(threshold=sys.maxsize) np.set_printoptions(linewidth=250) # example def main(): mesh_file = 'contact.msh' material_file = 'material.dat' aka.parseInput(material_file) spatial_dimension = 2 mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) # initialize model model = aka.SolidMechanicsModel(mesh) model.initFull(_analysis_method=aka._implicit_dynamic) # boundary conditions displacements = np.zeros(mesh.getNbNodes()*spatial_dimension) model.applyBC(aka.FixedValue(0., aka._x), 'lower_bottom') model.applyBC(aka.FixedValue(0., aka._y), 'lower_bottom') # Dirichlet - # nodes_top = mesh.getElementGroup('upper_top').getNodeGroup().getNodes().ravel() - # displacements = displacements.reshape([-1, 2]) - # displacements[nodes_top, 1] = -0.1 - # displacements = displacements.ravel() + model.applyBC(aka.FixedValue(-0.1, aka._y), 'upper_top') # to block the nodes + nodes_top = mesh.getElementGroup('upper_top').getNodeGroup().getNodes().ravel() + displacements = displacements.reshape([-1, 2]) + displacements[nodes_top, 1] = -0.1 + displacements = displacements.ravel() # Neumann (K is not invertible) - traction = np.zeros(spatial_dimension) - traction[1] = -1e9 - model.applyBC(aka.FromTraction(traction), 'upper_top') + # traction = np.zeros(spatial_dimension) + # traction[1] = -1e9 + # model.applyBC(aka.FromTraction(traction), 'upper_top') # init and solve model, data = init_model(model, mesh, mesh_file, material_file, spatial_dimension, displacements, 'lower_top', 'upper_bottom') - solve_step_direct(model, data, nb_max_iter=1, plot=True) + solve_step_direct(model, data, nb_max_iter=10, plot=True) def solve_step_direct(model, data, nb_max_iter=10, plot=False): #---------- assemble stiffness matrices: K, K_free ---------- model.assembleStiffnessMatrix() K_aka = model.dof_manager.getMatrix("K") K = sp.sparse.csc_matrix(aka.AkantuSparseMatrix(K_aka)) # K for all non blocked dofs rescaling = 30e9 # with E K_free = K[np.ix_(data['free_dofs'], data['free_dofs'])] / rescaling #---------- assemble external forces: f, f_ordering ---------- f = model.getExternalForce().ravel() # f for all non blocked dofs f_free = f[data['free_dofs']] #---------- initialize interface ---------- # 1i denotes interface 1, 2i denotes inteface 2 # take all boundary initialy as interface nodes # possible selection with radius calculation nodes1i = data['nodes1b'] positions1i = data['positions1b'] nodes2i = data['nodes2b'] positions2i = data['positions2b'] if plot: plt.figure() plt.triplot(data['positions'][:, 0], data['positions'][:, 1], data['connectivity']) plt.title('mesh') plt.xlabel('x') plt.ylabel('y') plt.axis('scaled') plt.title('Initial') # plt.ylim([0.9, 1.1]) plt.show() #---------- internodes iterations ---------- # run until converged or nb_max_iter is attained for i in range(nb_max_iter): print('--- iteration', i+1, '---') # select nodes belonging to interface nodes1i, nodes2i, positions1i, positions2i, radiuses1, radiuses2 = find_contact_nodes(nodes1i, nodes2i, positions1i, positions2i) dofs1i = nodes_to_dofs(nodes1i).ravel() dofs2i = nodes_to_dofs(nodes2i).ravel() nb_constraint_dofs = len(dofs1i) # global index of the interface among the boundary dofs sorter_free_dofs = np.argsort(data['free_dofs']) indx1i = sorter_free_dofs[np.searchsorted(data['free_dofs'], dofs1i, sorter=sorter_free_dofs)] indx2i = sorter_free_dofs[np.searchsorted(data['free_dofs'], dofs2i, sorter=sorter_free_dofs)] # subassemble matrices R12_normal, R21_normal, R12, R21 = assemble_Rijs(positions1i, positions2i, radiuses1, radiuses2) M1i, M2i = assemble_interface_masses(model, dofs1i, dofs2i) B, B_tilde, C = assemble_Bs(M1i, M2i, R12, R21, indx1i, indx2i, data['nb_free_dofs'], nb_constraint_dofs) A = assemble_A_explicit(K_free, B, B_tilde, C) b = assemble_b(f_free, R12_normal, K, data['displacements'], data['free_dofs'], data['blocked_dofs'], positions1i, positions2i, nb_constraint_dofs, rescaling) # solve positions_new, displacements, lambdas = solve_direct(A, b, data['positions'], data['displacements'], data['free_dofs'], nb_constraint_dofs, data['nb_dofs'], rescaling) if plot: plt.figure() plt.triplot(positions_new[:, 0], positions_new[:, 1], data['connectivity']) plt.title('mesh') plt.xlabel('x') plt.ylabel('y') plt.axis('scaled') plt.title('Iteration ' + str(i+1)) # plt.ylim([0.9, 1.1]) plt.show() # add or remove nodes nodes1i, nodes2i, diff_nb_nodes1i, diff_nb_nodes2i = remove_traction(positions_new, data['connectivity1b'], data['connectivity2b'], data['connectivity1b_body'], data['connectivity2b_body'], nodes1i, nodes2i, data['nodes1b'], data['nodes2b'], lambdas, R12, R21) positions1i = data['positions'][nodes1i, :] positions2i = data['positions'][nodes2i, :] if np.abs(diff_nb_nodes1i)+np.abs(diff_nb_nodes2i) == 0: print() print('successfully converged in', i+1, 'iterations') break return positions_new, displacements if __name__ == '__main__': main() diff --git a/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/example_iterative.py b/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/example_iterative.py index 9a45f3741..0817485b4 100644 --- a/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/example_iterative.py +++ b/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/example_iterative.py @@ -1,205 +1,207 @@ #!/usr/bin/env python # coding: utf-8 import akantu as aka import numpy as np import scipy as sp import matplotlib.pyplot as plt from scipy.sparse.linalg import LinearOperator from sksparse.cholmod import cholesky import sys sys.path.append("..") from prototype_internodes.functions import * from prototype_internodes.functions_contact_probl import * from prototype_internodes.init_model import init_model # example def main(): mesh_file = 'contact.msh' material_file = 'material.dat' aka.parseInput(material_file) spatial_dimension = 2 mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) # initialize model model = aka.SolidMechanicsModel(mesh) model.initFull(_analysis_method=aka._implicit_dynamic) # boundary conditions displacements = np.zeros(mesh.getNbNodes()*spatial_dimension) model.applyBC(aka.FixedValue(0., aka._x), 'lower_bottom') model.applyBC(aka.FixedValue(0., aka._y), 'lower_bottom') # Dirichlet - # nodes_top = mesh.getElementGroup('upper_top').getNodeGroup().getNodes().ravel() - # displacements = displacements.reshape([-1, 2]) - # displacements[nodes_top, 1] = -0.1 - # displacements = displacements.ravel() + model.applyBC(aka.FixedValue(-0.1, aka._y), 'upper_top') # to block the nodes + nodes_top = mesh.getElementGroup('upper_top').getNodeGroup().getNodes().ravel() + displacements = displacements.reshape([-1, 2]) + displacements[nodes_top, 1] = -0.1 + displacements = displacements.ravel() # Neumann (K is not invertible) - traction = np.zeros(spatial_dimension) - traction[1] = -1e9 - model.applyBC(aka.FromTraction(traction), "upper_top") + # traction = np.zeros(spatial_dimension) + # traction[1] = -1e9 + # model.applyBC(aka.FromTraction(traction), 'upper_top') # init and solve model, data = init_model(model, mesh, mesh_file, material_file, spatial_dimension, - displacements, "upper_bottom", "lower_top") + displacements, 'lower_top', 'upper_bottom') solve_step_iterative(model, data, nb_max_iter=10, plot=True) def solve_step_iterative(model, data, nb_max_iter=10, plot=False): #---------- assemble stiffness matrices: K, K_free ---------- model.assembleStiffnessMatrix() K_aka = model.dof_manager.getMatrix("K") K = sp.sparse.csc_matrix(aka.AkantuSparseMatrix(K_aka)) # K for all non blocked dofs rescaling = 30e9 # with E K_free = K[np.ix_(data['free_dofs'], data['free_dofs'])] / rescaling # calculate perturbed K free epsilon = 1e-8 K_tilde = K_free + epsilon*sp.sparse.identity(data['nb_free_dofs']) #---------- ordering and indexing of K and K_tilde ---------- sorter_free_dofs = np.argsort(data['free_dofs']) indx1b = sorter_free_dofs[np.searchsorted(data['free_dofs'], data['dofs1b'], sorter=sorter_free_dofs)] indx2b = sorter_free_dofs[np.searchsorted(data['free_dofs'], data['dofs2b'], sorter=sorter_free_dofs)] indxb = np.append(indx1b, indx2b) # ordering for cholesky of K_tilde with scikit-sparse factor = cholesky(K_tilde, beta=0) permut = factor.P() permut_inv = np.argsort(permut) indx_permut = np.setdiff1d(permut, indxb, assume_unique=True) ordering = np.append(indx_permut, indxb) ordering_inv = np.argsort(ordering) K_tilde_order = K_tilde[np.ix_(ordering, ordering)] K_free_order = K_free[np.ix_(ordering, ordering)] #---------- cholesky decomposition of K_tilde with scikit-sparse ---------- # this is not very efficient in sksparse chol_factor = cholesky(K_tilde_order, beta=0, ordering_method="natural") L = chol_factor.L() L = sp.sparse.csr_matrix(L) # explicit calculations of blocks boundary_start = data['nb_free_dofs'] - len(indxb) L22 = L[boundary_start:, boundary_start:] L22L22_T = L22.dot(L22.T) #---------- assemble external forces: f, f_ordering ---------- f = model.getExternalForce().ravel() # f for all non blocked dofs f_free = f[data['free_dofs']] f_free_order = f_free[ordering] #---------- initialize interface ---------- # 1i denotes interface 1, 2i denotes inteface 2 # take all boundary initialy as interface nodes # possible selection with radius calculation nodes1i = data['nodes1b'] positions1i = data['positions1b'] nodes2i = data['nodes2b'] positions2i = data['positions2b'] if plot: plt.figure() plt.triplot(data['positions'][:, 0], data['positions'][:, 1], data['connectivity']) plt.title('mesh') plt.xlabel('x') plt.ylabel('y') plt.axis('scaled') plt.title('Initial') # plt.ylim([0.9, 1.1]) plt.show() #---------- internodes iterations ---------- # run until converged or nb_max_iter is attained for i in range(nb_max_iter): print('--- iteration', i+1, '---') # select nodes belonging to interface nodes1i, nodes2i, positions1i, positions2i, radiuses1, radiuses2 = find_contact_nodes(nodes1i, nodes2i, positions1i, positions2i) dofs1i = nodes_to_dofs(nodes1i).ravel() dofs2i = nodes_to_dofs(nodes2i).ravel() nb_constraint_dofs = len(dofs1i) # local index of the interface among the boundary dofs sorter_boundary_dofs = np.argsort(data['boundary_dofs']) indx1i_loc = sorter_boundary_dofs[np.searchsorted(data['boundary_dofs'], dofs1i, sorter=sorter_boundary_dofs)] indx2i_loc = sorter_boundary_dofs[np.searchsorted(data['boundary_dofs'], dofs2i, sorter=sorter_boundary_dofs)] # global index of the interface among all dofs boundary_start = data['nb_free_dofs'] - len(indxb) indx1i = indx1i_loc + boundary_start indx2i = indx2i_loc + boundary_start # subassemble matrices R12_normal, R21_normal, R12, R21 = assemble_Rijs(positions1i, positions2i, radiuses1, radiuses2) M1i, M2i = assemble_interface_masses(model, dofs1i, dofs2i) B, B_tilde, C = assemble_Bs(M1i, M2i, R12, R21, indx1i, indx2i, data['nb_free_dofs'], nb_constraint_dofs) b = assemble_b(f_free_order, R12_normal, K, data['displacements'], data['free_dofs'], data['blocked_dofs'], positions1i, positions2i, nb_constraint_dofs, rescaling) + b[:data['nb_free_dofs']] = b[:data['nb_free_dofs']][ordering] A_size = data['nb_free_dofs'] + nb_constraint_dofs A_op = lambda x: linearoperator_A(x, K_tilde_order, B, B_tilde, data['nb_free_dofs']) A_op_spla = LinearOperator((A_size, A_size), matvec=A_op) # prepare preconditioner V = prepare_precond(M1i, M2i, R12, R21, L22, L22L22_T, indx1i_loc, indx2i_loc, data['nb_boundary_dofs']) M_inv_op = lambda y: linearoperator_precond(y, L22, V, B, M1i, chol_factor, data['nb_free_dofs'], data['nb_boundary_dofs']) M_inv_op_spla = LinearOperator((A_size, A_size), matvec=M_inv_op) # solve positions_new, displacements, lambdas = solve_iterative(A_op_spla, b, M_inv_op_spla, data['positions'], data['displacements'], data['free_dofs'], ordering, nb_constraint_dofs, data['nb_dofs'], rescaling) if plot: plt.figure() plt.triplot(positions_new[:, 0], positions_new[:, 1], data['connectivity']) plt.title('mesh') plt.xlabel('x') plt.ylabel('y') plt.axis('scaled') plt.title('Iteration ' + str(i+1)) # plt.ylim([0.9, 1.1]) plt.show() # add or remove nodes nodes1i, nodes2i, diff_nb_nodes1i, diff_nb_nodes2i = remove_traction(positions_new, data['connectivity1b'], data['connectivity2b'], data['connectivity1b_body'], data['connectivity2b_body'], nodes1i, nodes2i, data['nodes1b'], data['nodes2b'], lambdas, R12, R21) positions1i = data['positions'][nodes1i, :] positions2i = data['positions'][nodes2i, :] if np.abs(diff_nb_nodes1i)+np.abs(diff_nb_nodes2i) == 0: print() print('successfully converged in', i+1, 'iterations') break return positions_new, displacements if __name__ == '__main__': main() diff --git a/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/functions.py b/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/functions.py index 2df9ec4ad..b1b6794b9 100644 --- a/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/functions.py +++ b/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/functions.py @@ -1,288 +1,285 @@ #!/usr/bin/env python # coding: utf-8 import akantu as aka import numpy as np import scipy as sp import scipy.sparse.linalg as spla from scipy.sparse.linalg import spsolve def nodes_to_dofs(nodes): """Returns dof numbers for nodes in 2D. """ return np.column_stack([2*nodes, 2*nodes+1]).reshape([-1, 1]) def assemble_interface_masses(model, dofs1i, dofs2i): """Assembles mass matrices for interfaces. :param model: akantu solid mechanics model :param dofs1i: dofs of body 1 interface (master) :param dofs2i: dofs of body 2 interface (slave) :returns: sparse matrices of interface mass """ # assemble and store mass model.assembleMass() M_aka = model.getDOFManager().getMatrix('M') M = sp.sparse.lil_matrix(aka.AkantuSparseMatrix(M_aka)) # select only the dofs of the interface M1i = M[np.ix_(dofs1i, dofs1i)] M2i = M[np.ix_(dofs2i, dofs2i)] M1i = sp.sparse.csc_matrix(M1i) M2i = sp.sparse.csc_matrix(M2i) return M1i, M2i def assemble_Bs(M1i, M2i, R12, R21, indx1i, indx2i, nb_free_dofs, nb_constraint_dofs): """Assembles B and B_tilde for internodes matrix. :param M1i: sparse matrix interface mass 1 :param M2i: sparse matrix interface mass 2 :param R12: interpolation matrix master to slave :param R21: interpolation matrix slave to master :param indx1i: global indices of interface 1 (master) :param indx2i: global indices of interface 2 (slave) :param nb_free_dofs: number of non blocked dofs :param nb_constraint_dofs: number of master dofs which act as constraint :returns: sparse blocke matrices for internodes formulation """ - B = sp.sparse.lil_matrix(np.zeros((nb_free_dofs, nb_constraint_dofs))) B_tilde = sp.sparse.lil_matrix(np.zeros((nb_constraint_dofs, nb_free_dofs))) C = sp.sparse.lil_matrix(np.zeros((nb_constraint_dofs, nb_constraint_dofs))) B[indx1i, :] = - M1i B[indx2i, :] = M2i * R21 B_tilde[:, indx1i] = sp.sparse.eye(nb_constraint_dofs) B_tilde[:, indx2i] = - R12 return B, B_tilde, C def assemble_A_explicit(K_free, B, B_tilde, C): """Explicitly assembles internodes matrix. :param K_free: stiffness matrix with non blocked dofs :param B: B matrix of internodes formulation :param B_tilde: B_tilde matrix of internodes formulation :param C: all zero matrix of internodes formulation """ A1 = sp.sparse.hstack([K_free, B]) A2 = sp.sparse.hstack([B_tilde, C]) A = sp.sparse.vstack([A1, A2]) A = sp.sparse.csr_matrix(A) return A def assemble_b(f_free, R12_normal, K, displacements, free_dofs, blocked_dofs, positions1i, positions2i, nb_constraint_dofs, rescaling): """Assemble right hand side of Ax=b. :param f_free: force vector of free dofs :param R12_normal: nodal interpolation matrix master to slave :param K: stiffness matrix + :param displacements: displacements vector :param free_dofs: dofs numbers of non blocked dofs :param blocked_dofs: dofs numbers of blocked dofs - :param positions1i: initial positions of interface 1 nodes - :param positions2i: initial positions of interface 2 nodes + :param positions1i: initial positions of interface 1 nodes + :param positions2i: initial positions of interface 2 nodes :param nb_constraint_dofs: number of master dofs which act as constraint - :param nb_dofs: total number of dofs :param rescaling: rescaling stiffness matrix """ nb_free_dofs = len(free_dofs) b = np.zeros(nb_free_dofs + nb_constraint_dofs) # Dirichlet displacements K_blocked = K[free_dofs, :] K_blocked = K_blocked[:, blocked_dofs] dirichlet = K_blocked.dot(displacements[blocked_dofs]) b[0:nb_free_dofs] = 1/rescaling * (f_free - dirichlet) b[nb_free_dofs:] = (R12_normal.dot(positions2i) - positions1i).ravel() return b def solve_direct(A, b, positions, displacements, free_dofs, nb_constraint_dofs, nb_dofs, rescaling): """Direct solve of internodes equation Ax=b. :param A: sparse matrix internodes :param b: right hand side vector - :param displacements: displacements vector :param positions: initial positions of all nodes - :param scaling: scaling factor for stiffness matrix + :param displacements: displacements vector :param free_dofs: dofs numbers of non blocked dofs :param nb_constraint_dofs: number of master dofs which act as constraint :param nb_dofs: total number of dofs :param rescaling: rescaling stiffness matrix :returns: new positions, displacements and lambdas resulting from constraints """ nb_free_dofs = len(free_dofs) x = spsolve(A, b) displacements[free_dofs] = x[:nb_free_dofs] positions_new = positions + displacements.reshape([-1, 2]) lambdas = rescaling * x[nb_free_dofs:nb_free_dofs+nb_constraint_dofs].reshape([-1, 2]) return positions_new, displacements, lambdas # all functions below are needed for iterative solve def solve_iterative(A_op_spla, b, M_inv_op_spla, positions, displacements, free_dofs, ordering, nb_constraint_dofs, nb_dofs, rescaling): """Iterative solve of internodes equation Ax=b. :param A_op_spla: scipy linear operator for A :param b: right hand side vector - :param M_inv_op_spla: scipy linear operator for B + :param M_inv_op_spla: scipy linear operator for M :param positions: initial positions of all nodes + :param displacements: displacements vector + :param free_dofs: dofs numbers of non blocked dofs :param ordering: reordered stiffness matrix due to cholesky factorization - :param scaling: scaling factor for stiffness matrix - :param free_dofs: mask of free dofs :param nb_constraint_dofs: number of master dofs which act as constraint :param nb_dofs: total number of dofs + :param rescaling: rescaling stiffness matrix :returns: new positions and lambdas resulting from constraints """ nb_free_dofs = len(free_dofs) x, exitCode = spla.gmres(A_op_spla, b, tol=1e-07, restart=5, maxiter=3, M=M_inv_op_spla) disp_order = x[:nb_free_dofs] inv_ordering = np.argsort(ordering) displacements[free_dofs] = disp_order[inv_ordering] positions_new = positions + displacements.reshape([-1, 2]) lambdas = rescaling * x[nb_free_dofs:nb_free_dofs+nb_constraint_dofs].reshape([-1, 2]) return positions_new, displacements, lambdas def linearoperator_A(x, K_free_order, B, B_tilde, nb_free_dofs): """Creates linear operator for internodes matrix. - :param K_free: stiffness matrix with non blocked dofs + :param K_free_order: ordered stiffness matrix with non blocked dofs :param B: B matrix of internodes formulation :param B_tilde: B_tilde matrix of internodes formulation :param nb_free_dofs: number of non blocked dofs - :param nb_free_dofs: number of non blocked dofs :param nb_constraint_dofs: number of master dofs which act as constraint """ # A @ [x1; x2] = [y1; y2] x1 = x[:nb_free_dofs] x2 = x[nb_free_dofs:] y1 = K_free_order.dot(x1) + B.dot(x2) y2 = B_tilde.dot(x1) y = np.concatenate((y1, y2)) return y def prepare_precond(M1i, M2i, R12, R21, L22, L22L22_T, indx1i, indx2i, nb_boundary_dofs): """Prepares matrix V needed for construction of preconditioner. :param M1i: sparse matrix interface mass 1 :param M2i: sparse matrix interface mass 2 :param R12: interpolation matrix master to slave :param R21: interpolation matrix slave to master :param L22: block 22 of cholesky factorization :param L22L22_T: precomputed L_22 * L_22^T :param indx1i: indices of interface 1 (master) :param indx2i: indices of interface 2 (slave) - :param nb_boundary_dofs: indices of initial boundary interface dofs + :param nb_boundary_dofs: number of initial boundary dofs :return: matrix V needed for preconditioner """ - indxi = np.append(indx1i, indx2i) # assemble Y1*Y2^T Q1 = -M2i.dot((R21.dot(spla.inv(M1i)))) Q2 = -R12.T Y1 = sp.sparse.vstack([sp.sparse.identity(len(indx1i)), Q1]) Y2 = sp.sparse.vstack([sp.sparse.identity(len(indx1i)), Q2]) Y1Y2_T = Y1.dot(Y2.T) Z = sp.sparse.lil_matrix((nb_boundary_dofs, nb_boundary_dofs)) Z[np.ix_(indxi, indxi)] = Y1Y2_T # assemble V = L22*L22^T - alpha*Y1*Y2^T alpha = 7.922 # fixed for preconditioner V = L22L22_T - alpha*Z return V def linearoperator_precond(y, L22, V, B, M1i, chol_factor, nb_free_dofs, nb_boundary_dofs): """Creates a linear operator for inverse of preconditioner M^-1 , where MAx=Mb and Ax=y. :param y: solution to solve for :param L22: submatrice of chelosky for perturbed stiffness matrix :param V: prepared matrix for preconditioner :param B: B matrix of internodes formulation - :param M1i: sparse matrix interface mass 1 + :param M1i: sparse matrix interface mass 1 :param chol_factor: sksparse factor to do operations with L :param nb_free_dofs: number of non blocked dofs :param nb_boundary_dofs: indices of initial boundary interface dofs """ - boundary_start = nb_free_dofs - nb_boundary_dofs # M @ [x1; x2] = [y1; y2] y1 = y[:nb_free_dofs] y2 = y[nb_free_dofs:] alpha = 7.922 # fixed for preconditioner x2 = spla.spsolve(M1i, -alpha*y2) y1_tilde = y1 - 2*B.dot(x2) x1 = chol_factor.solve_L(y1_tilde, use_LDLt_decomposition=False) x1a = x1[:boundary_start] x1b = L22.dot(x1[boundary_start:]) # solving for V could be done # with GMRES as well + preconditiong x1b = spla.spsolve(V, x1b) x1b = L22.T.dot(x1b) x1 = np.concatenate((x1a, x1b)) x1 = chol_factor.solve_Lt(x1, use_LDLt_decomposition=False) x = np.concatenate((x1, x2)) return x def solve_iterative(A_op_spla, b, M_inv_op_spla, positions, displacements, free_dofs, ordering, nb_constraint_dofs, nb_dofs, rescaling): """Iterative solve of internodes equation Ax=b. :param A_op_spla: scipy linear operator for A :param b: right hand side vector :param M_inv_op_spla: scipy linear operator for B :param positions: initial positions of all nodes + :param displacements: displacements vector + :param free_dofs: dofs numbers of non blocked dofs :param ordering: reordered stiffness matrix due to cholesky factorization - :param scaling: scaling factor for stiffness matrix - :param free_dofs: mask of free dofs :param nb_constraint_dofs: number of master dofs which act as constraint :param nb_dofs: total number of dofs + :param rescaling: rescaling stiffness matrix :returns: new positions and lambdas resulting from constraints """ nb_free_dofs = len(free_dofs) x, exitCode = spla.gmres(A_op_spla, b, tol=1e-07, restart=5, maxiter=3, M=M_inv_op_spla) disp_order = x[:nb_free_dofs] inv_ordering = np.argsort(ordering) displacements[free_dofs] = disp_order[inv_ordering] positions_new = positions + displacements.reshape([-1, 2]) lambdas = rescaling * x[nb_free_dofs:nb_free_dofs+nb_constraint_dofs].reshape([-1, 2]) return positions_new, displacements, lambdas diff --git a/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/functions_contact_probl.py b/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/functions_contact_probl.py index ed68b58fa..2568f420d 100644 --- a/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/functions_contact_probl.py +++ b/test/test_model/test_contact_mechanics_internodes_model/prototype_internodes/functions_contact_probl.py @@ -1,252 +1,321 @@ #!/usr/bin/env python # coding: utf-8 import numpy as np import scipy as sp from scipy.linalg import lstsq from scipy import spatial import sys sys.path.append("..") from prototype_internodes.functions import nodes_to_dofs -def find_contact_nodes(nodes1i, nodes2i, coords1i, coords2i): - radiuses1, nnzR21 = compute_radiuses(nodes1i, nodes2i, coords1i, coords2i) - radiuses2, nnzR12 = compute_radiuses(nodes2i, nodes1i, coords2i, coords1i) +def find_contact_nodes(nodes1i, nodes2i, positions1i, positions2i): + """Contact nodes algorithm. + + :param nodes1i: nodes of body 1 interface (master) + :param nodes2i: nodes of body 2 interface (slave) + :param positions1i: initial positions of body 1 interface (master) + :param positions2i: initial positions of body 2 interface (slave) + :returns: selected nodes and radiuses + """ + radiuses1, nnzR21 = compute_radiuses(nodes1i, nodes2i, positions1i, positions2i) + radiuses2, nnzR12 = compute_radiuses(nodes2i, nodes1i, positions2i, positions1i) nodes1i_mask = nnzR12 > 0 nodes2i_mask = nnzR21 > 0 while np.any(nodes1i_mask == False) or np.any(nodes2i_mask == False): nodes1i = nodes1i[nodes1i_mask] nodes2i = nodes2i[nodes2i_mask] - coords1i = coords1i[nodes1i_mask] - coords2i = coords2i[nodes2i_mask] + positions1i = positions1i[nodes1i_mask] + positions2i = positions2i[nodes2i_mask] dofs1i = nodes_to_dofs(nodes1i).ravel() dofs2i = nodes_to_dofs(nodes2i).ravel() - radiuses1, nnzR21 = compute_radiuses(nodes1i, nodes2i, coords1i, coords2i) - radiuses2, nnzR12 = compute_radiuses(nodes2i, nodes1i, coords2i, coords1i) + radiuses1, nnzR21 = compute_radiuses(nodes1i, nodes2i, positions1i, positions2i) + radiuses2, nnzR12 = compute_radiuses(nodes2i, nodes1i, positions2i, positions1i) nodes1i_mask = nnzR12 > 0 nodes2i_mask = nnzR21 > 0 - return nodes1i, nodes2i, coords1i, coords2i, radiuses1, radiuses2 + return nodes1i, nodes2i, positions1i, positions2i, radiuses1, radiuses2 -def compute_radiuses(nodes1i, nodes2i, coords1i, coords2i): +def compute_radiuses(nodes1i, nodes2i, positions1i, positions2i): + """Computes the radiuses of attack. + + :param nodes1i: nodes of body 1 interface (master) + :param nodes2i: nodes of body 2 interface (slave) + :param positions1i: initial positions of body 1 interface (master) + :param positions2i: initial positions of body 2 interface (slave) + :returns: radiuses, and resulting nonzero components of R matrix + """ + # fixed parameters c = 0.5 # conditition (2) C = 0.95 # condition (3) - n = 1 # consider n nearest neighboors + n = 1 # consider 1 nearest neighboors d = 0.05 # tolerance, for radius of "attack" estimation - M = len(coords1i) - N = len(coords2i) + M = len(positions1i) + N = len(positions2i) radiuses = np.zeros(M) nnzRMM = np.zeros(M) nnzRNM = np.zeros(N) nnzCMM = np.zeros(M) nnzCNM = np.zeros(M) maxS = np.inf f = 0 niter = 0 maxiter = 10 while maxS > f and niter < maxiter-1: f = np.floor(1/(np.power(1-c, 4)*(1+4*c))) # maximum number of supports for k in range(M): - point = coords1i[k, :].reshape(1, -1) - neighbors = coords1i.copy() + point = positions1i[k, :].reshape(1, -1) + neighbors = positions1i.copy() neighbors[k, :] = np.inf distMM = spatial.distance.cdist(neighbors, point).ravel() - distMN = spatial.distance.cdist(coords2i, point).ravel() + distMN = spatial.distance.cdist(positions2i, point).ravel() rMM = np.min(distMM) rNM = np.sqrt(d*d + 0.25*np.power(rMM, 2)) radius = np.maximum(rMM, rNM) - # rMM = distMM[np.argpartition(distMM, n)[:n]] - # rNM = np.sqrt(d*d + 0.25*np.power(rMM, 2)) - # radius = np.maximum([rMM[-1], rNM]) - - # if radius > rMM[0]/c: - # radius = rMM[0]/c - if radius > rMM/c: radius = rMM/c s1 = distMM < radius s2 = distMN < C*radius nnzRMM[s1] = nnzRMM[s1] + 1 nnzRNM[s2] = nnzRNM[s2] + 1 nnzCMM[k] = np.sum(s1) nnzCNM[k] = np.sum(s2) radiuses[k] = radius maxS = np.max(nnzRMM) if maxS > f: c = 0.5 * (1 + c); nnzRMM = np.zeros(M) nnzRNM = np.zeros(N) nnzCMM = np.zeros(M) nnzCNM = np.zeros(M) niter = niter+1 return radiuses, nnzRNM def wendland(dists, radiuses): + """ Compute the Beckert & Wendland RBF + + :param dists: distances + :param radiuses: radius for each distance (same length) + """ result = np.zeros(len(dists)) mask = dists <= radiuses result[mask] = np.power(1-dists[mask]/radiuses[mask], 4) * (1+4*dists[mask]/radiuses[mask]) return result -def phi_constructor(coords_i, coords_j, radiuses_j, rad_func): - N = len(coords_i) - M = len(coords_j) +def phi_constructor(positions_i, positions_j, radiuses_j, rad_func): + """Construct Phi for RBF interpolation. + + :param positions_i: positions of evaluation points + :param positions_j: positions of reference points + :param radiuses_j: radiuses of reference points + """ + N = len(positions_i) + M = len(positions_j) - dists = spatial.distance.cdist(coords_i, coords_j) + dists = spatial.distance.cdist(positions_i, positions_j) radiuses_j = np.tile(radiuses_j, N) phi = rad_func(dists.ravel(), radiuses_j.ravel()) return phi.reshape([N, M]) -def Rij_constructor(coords_i, coords_j, radiuses_j): - phiMM = phi_constructor(coords_j, coords_j, radiuses_j, wendland) - phiNM = phi_constructor(coords_i, coords_j, radiuses_j, wendland) +def Rij_constructor(positions_i, positions_j, radiuses_j): + """Construct a RBF interpolation. + + :param positions_i: positions of evaluation points + :param positions_j: positions of reference points + :param radiuses_j: radiuses of reference points + """ + phiMM = phi_constructor(positions_j, positions_j, radiuses_j, wendland) + phiNM = phi_constructor(positions_i, positions_j, radiuses_j, wendland) Rij = phiNM.dot(np.linalg.inv(phiMM)) g = Rij.dot(np.ones((Rij.shape[1], 1))) Rij_norm = Rij * (1/g) return Rij_norm -def assemble_Rijs(coords1i, coords2i, radiuses1, radiuses2): - R12_normal = Rij_constructor(coords1i, coords2i, radiuses2) - R21_normal = Rij_constructor(coords2i, coords1i, radiuses1) +def assemble_Rijs(positions1i, positions2i, radiuses1, radiuses2): + """ Assembles the Radial Basis function (RBF) interpolation matrices. + + :param positions_i: positions of evaluation points + :param positions_j: positions of reference points + :param radiuses_i: radiuses of evaluation points + :param radiuses_j: radiuses of reference points + """ + R12_normal = Rij_constructor(positions1i, positions2i, radiuses2) + R21_normal = Rij_constructor(positions2i, positions1i, radiuses1) R21 = sp.sparse.csr_matrix(extend_to_2D(R21_normal)) R12 = sp.sparse.csr_matrix(extend_to_2D(R12_normal)) return R12_normal, R21_normal, R12, R21 def extend_to_2D(R): + """ Extend the RBF matrices to 2D. """ R_extended = np.repeat(np.repeat(R,2,axis=1), 2, axis=0) R_extended[1::2,::2] = 0 R_extended[::2,1::2] = 0 return R_extended -def remove_traction(coords_new, connectivity1b, connectivity2b, connectivity1b_body, connectivity2b_body, +def remove_traction(positions_new, connectivity1b, connectivity2b, connectivity1b_body, connectivity2b_body, nodes1i, nodes2i, nodes1b, nodes2b, lambda1, R12, R21): - normals1b = compute_normals(coords_new, nodes1b, connectivity1b, connectivity1b_body) - normals2b = compute_normals(coords_new, nodes2b, connectivity2b, connectivity2b_body) + + """ Check if there is still traction between the bodies. + + :param positions_new: new positions after internodes step + :param connectivity1b: connectivity of boundary surface 1 (master) + :param connectivity2b: connectivity of boundary surface 2 (slave) + :param connectivity1b_body: connectivity of boundary elements 1 (master) + :param connectivity2b_body: connectivity of boundary elements 2 (slave) + :param nodes1i: nodes of interface 1 + :param nodes2i: nodes of interface 2 + :param nodes1b: nodes of boundary 1 (potential interface nodes) + :param nodes2b: nodes of boundary 2 (potential interface nodes) + :param lambda1: lambdas of interface 1 (solution) + :param R12: RBF interpolation master to slave + :param R21: RBF interpolation slave to master + """ + normals1b = compute_normals(positions_new, nodes1b, connectivity1b, connectivity1b_body) + normals2b = compute_normals(positions_new, nodes2b, connectivity2b, connectivity2b_body) normals1i = normals1b[np.in1d(nodes1b, nodes1i)] normals2i = normals2b[np.in1d(nodes2b, nodes2i)] lambda2 = -R21.dot(lambda1.reshape([-1, 1])).reshape([-1, 2]) scalar1 = np.sum(lambda1*normals1i, axis=1) scalar2 = np.sum(lambda2*normals2i, axis=1) nodes1i_dump = nodes1i[scalar1>0] nodes2i_dump = nodes2i[scalar2>0] if len(nodes1i_dump) == 0 and len(nodes2i_dump) == 0: # gap verification - nodes1i_add, nodes2i_add = detect_gaps(coords_new, nodes1i, nodes2i, normals1i, normals2i) + nodes1i_add, nodes2i_add = detect_gaps(positions_new, nodes1i, nodes2i, normals1i, normals2i) nodes1i, diff_nb_nodes1i = update_interface(nodes1i_add, nodes1i, 'add') nodes2i, diff_nb_nodes2i = update_interface(nodes2i_add, nodes2i, 'add') + + print(diff_nb_nodes1i, ' nodes added to interface 1') + print(diff_nb_nodes2i, ' nodes added to interface 2') else: nodes1i, diff_nb_nodes1i = update_interface(nodes1i_dump, nodes1i, 'dump') nodes2i, diff_nb_nodes2i = update_interface(nodes2i_dump, nodes2i, 'dump') - print(diff_nb_nodes1i, ' nodes removed from interface 1') - print(diff_nb_nodes2i, ' nodes removed from interface 2') + print(diff_nb_nodes1i, ' nodes removed from interface 1') + print(diff_nb_nodes2i, ' nodes removed from interface 2') return nodes1i, nodes2i, diff_nb_nodes1i, diff_nb_nodes2i def update_interface(new_nodes, nodesi, case): + """ Remove of add interfaces after traction check. """ if case == 'dump': nodesi_new = nodesi[~np.in1d(nodesi, new_nodes)] if case == 'add': nodesi_new = np.union1d(nodesi, new_nodes) diff_nb_nodes = len(nodesi) -len(nodesi_new) return nodesi_new, diff_nb_nodes -def compute_normals(coords_new, nodesb, connectivityb, connectivityb_body): +def compute_normals(positions_new, nodesb, connectivityb, connectivityb_body): + """Comput normals on interface surface. + + :param positions_new: new positions after internodes step + :param nodesb: nodes of boundary(potential interface nodes) + :param connectivityb: connectivity of boundary surface + """ n = len(nodesb) m = len(connectivityb) connectivityi_body = connectivityb_body[np.in1d(connectivityb_body, nodesb).reshape(connectivityb_body.shape).any(axis=1)] nodesb_body = np.unique(connectivityi_body[~np.isin(connectivityb_body, nodesb)]) - tangents = coords_new[connectivityb[:, 1]] - coords_new[connectivityb[:, 0]] + tangents = positions_new[connectivityb[:, 1]] - positions_new[connectivityb[:, 0]] lengths = np.linalg.norm(tangents, axis=1).reshape([-1,1]) tangents = tangents/lengths normals = np.zeros((m, 2)) normals[:, 0] = -tangents[:, 1] normals[:, 1] = tangents[:, 0] normals_avg = np.zeros((n, 2)) gamma = 1e-3 # step size for j in range(n): node = nodesb[j] - coord = coords_new[node, :] + coord = positions_new[node, :] id = np.in1d(connectivityb, node).reshape(connectivityb.shape).any(axis=1) length = lengths[id] normal_avg = 1/np.sum(length)*np.sum(normals[id, :]*length, axis=0) tang_plus = (coord + gamma*normal_avg).reshape([-1, 2]) tang_minus = (coord - gamma*normal_avg).reshape([-1, 2]) - min_plus = np.min(spatial.distance.cdist(coords_new[nodesb_body, :], tang_plus).ravel()) - min_minus = np.min(spatial.distance.cdist(coords_new[nodesb_body, :], tang_minus).ravel()) + min_plus = np.min(spatial.distance.cdist(positions_new[nodesb_body, :], tang_plus).ravel()) + min_minus = np.min(spatial.distance.cdist(positions_new[nodesb_body, :], tang_minus).ravel()) if min_plus > min_minus: normals_avg[j, :] = normal_avg else: normals_avg[j, :] = -normal_avg norms = np.linalg.norm(normals_avg, axis=1).reshape([-1, 1]) normals_avg = (normals_avg/norms).reshape([-1, 2]) return normals_avg -def detect_gaps(coords_new, nodes1i, nodes2i, normals1i, normals2i): +def detect_gaps(positions_new, nodes1i, nodes2i, normals1i, normals2i): + """Detect gaps between interfaces. + + :param positions_new: new positions after internodes step + :param nodes1i: nodes of interface 1 + :param nodes2i: nodes of interface 2 + :param normals1i: normal vectors of nodes on interface 1 + :param normals2i: normal vectors of nodes on interface 2 + """ tol = 0.9 # tolerance for gap detection (could be changed as input) h = 0.05 # mesh size (shouldn't be fixed!) - coords1i = coords_new[nodes1i, :] - coords2i = coords_new[nodes2i, :] + positions1i = positions_new[nodes1i, :] + positions2i = positions_new[nodes2i, :] - nodes1i, nodes2i, coords1i, coords2i, radiuses1, radiuses2 = find_contact_nodes(nodes1i, nodes2i, coords1i, coords2i) + nodes1i, nodes2i, positions1i, positions2i, radiuses1, radiuses2 = find_contact_nodes(nodes1i, nodes2i, positions1i, positions2i) - R21_normal = Rij_constructor(coords2i, coords1i, radiuses1) + R21_normal = Rij_constructor(positions2i, positions1i, radiuses1) R21 = sp.sparse.csr_matrix(extend_to_2D(R21_normal)) - R12_normal = Rij_constructor(coords1i, coords2i, radiuses2) + R12_normal = Rij_constructor(positions1i, positions2i, radiuses2) R12 = sp.sparse.csr_matrix(extend_to_2D(R12_normal)) - diffs1 = R12.dot(coords2i.reshape([-1, 1])) - coords1i.reshape([-1, 1]) + diffs1 = R12.dot(positions2i.reshape([-1, 1])) - positions1i.reshape([-1, 1]) diffs1 = diffs1.reshape([-1, 2]) - diffs2 = R21.dot(coords1i.reshape([-1, 1])) - coords2i.reshape([-1, 1]) + diffs2 = R21.dot(positions1i.reshape([-1, 1])) - positions2i.reshape([-1, 1]) diffs2 = diffs2.reshape([-1, 2]) scalar1 = np.sum(diffs1*normals1i, axis=1) scalar2 = np.sum(diffs2*normals2i, axis=1) threshold = -tol*h nodes1i_add = nodes1i[scalar1