diff --git a/common/aka_math_inline_impl.cc b/common/aka_math_inline_impl.cc index ae01931d5..2b21051c9 100644 --- a/common/aka_math_inline_impl.cc +++ b/common/aka_math_inline_impl.cc @@ -1,373 +1,373 @@ /** * @file aka_math_inline_impl.cc * @author Nicolas Richart * @date Wed Jul 28 13:20:35 2010 * * @brief Implementation of the inline functions of the math toolkit * * @section LICENSE * * \ * */ #ifdef AKANTU_USE_CBLAS # ifndef AKANTU_USE_CBLAS_MKL # include # else // AKANTU_USE_CBLAS_MKL # include # endif //AKANTU_USE_CBLAS_MKL #endif //AKANTU_USE_CBLAS /* -------------------------------------------------------------------------- */ inline void Math::matrix_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y) { #ifdef AKANTU_USE_CBLAS /// y = alpha*op(A)*x + beta*y cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 1, A, n, x, 1, 0, y, 1); #else memset(y, 0, m*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * n; for (UInt j = 0; j < n; ++j) { y[i] += A[A_i + j] * x[j]; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1, A, k, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[A_i + l] * B[l * n + j]; } } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, k, 1, A, m, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { // UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[l * n + j]; } } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, n, k, 1, A, k, B, k, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[A_i + l] * B[B_j + l]; } } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, m, n, k, 1, A, m, B, k, 0, C, n); #else memset(C, 0, m * n * sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[B_j + l]; } } } #endif } /* -------------------------------------------------------------------------- */ inline Real Math::det2(const Real * mat) { return mat[0]*mat[3] - mat[1]*mat[2]; } /* -------------------------------------------------------------------------- */ inline Real Math::det3(const Real * mat) { return mat[0]*(mat[4]*mat[8]-mat[7]*mat[5]) - mat[3]*(mat[1]*mat[8]-mat[7]*mat[2]) + mat[6]*(mat[1]*mat[5]-mat[4]*mat[2]); } /* -------------------------------------------------------------------------- */ inline void Math::normal2(const Real * vec,Real * normal) { normal[0] = vec[1]; normal[1] = -vec[0]; Math::normalize2(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normal3(const Real * vec1,const Real * vec2,Real * normal) { Math::vectorProduct3(vec1,vec2,normal); Math::normalize3(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normalize2(Real * vec) { Real norm = Math::norm2(vec); vec[0] /= norm; vec[1] /= norm; } /* -------------------------------------------------------------------------- */ inline void Math::normalize3(Real * vec) { Real norm = Math::norm3(vec); vec[0] /= norm; vec[1] /= norm; vec[2] /= norm; } /* -------------------------------------------------------------------------- */ inline Real Math::norm2(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::norm3(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]); } /* -------------------------------------------------------------------------- */ inline void Math::inv2(const Real * mat,Real * inv) { Real det_mat = det2(mat); inv[0] = mat[3] / det_mat; inv[1] = -mat[1] / det_mat; inv[2] = -mat[2] / det_mat; inv[3] = mat[0] / det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::inv3(const Real * mat,Real * inv) { Real det_mat = det3(mat); inv[0] = (mat[4]*mat[8] - mat[7]*mat[5])/det_mat; inv[1] = (mat[2]*mat[7] - mat[8]*mat[1])/det_mat; inv[2] = (mat[1]*mat[5] - mat[4]*mat[2])/det_mat; inv[3] = (mat[5]*mat[6] - mat[8]*mat[3])/det_mat; inv[4] = (mat[0]*mat[8] - mat[6]*mat[2])/det_mat; inv[5] = (mat[2]*mat[3] - mat[5]*mat[0])/det_mat; inv[6] = (mat[3]*mat[7] - mat[6]*mat[4])/det_mat; inv[7] = (mat[1]*mat[6] - mat[7]*mat[0])/det_mat; inv[8] = (mat[0]*mat[4] - mat[3]*mat[1])/det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::vectorProduct3(const Real * v1, const Real * v2, Real * res) { res[0] = v1[1]*v2[2] - v1[2]*v2[1]; res[1] = v1[2]*v2[0] - v1[0]*v2[2]; res[2] = v1[0]*v2[1] - v1[1]*v2[0]; } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot2(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot3(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_2d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1])); } /* -------------------------------------------------------------------------- */ inline Real Math::triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3) { /** * @f{eqnarray*}{ * r &=& A / s \\ * A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)} \\ * s &=& \frac{a + b + c}{2} * @f} */ Real a, b, c; a = distance_2d(coord1, coord2); b = distance_2d(coord2, coord3); c = distance_2d(coord1, coord3); Real s; s = (a + b + c) * 0.5; return sqrt((s - a) * (s - b) * (s - c) / s); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_3d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1]) + (y[2] - x[2])*(y[2] - x[2]) ); } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real xx[9], vol; xx[0] = coord2[0]; xx[1] = coord2[1]; xx[2] = coord2[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol = det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol -= det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol += det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord3[0]; xx[7] = coord3[1]; xx[8] = coord3[2]; vol -= det3(xx); vol /= 6; return vol; } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real l12, l13, l14, l23, l24, l34; l12 = distance_3d(coord1, coord2); l13 = distance_3d(coord1, coord3); l14 = distance_3d(coord1, coord4); l23 = distance_3d(coord2, coord3); l24 = distance_3d(coord2, coord4); l34 = distance_3d(coord3, coord4); Real s1, s2, s3, s4; s1 = (l12 + l23 + l13) * 0.5; s1 = sqrt(s1*(s1-l12)*(s1-l23)*(s1-l13)); s2 = (l12 + l24 + l14) * 0.5; s2 = sqrt(s2*(s2-l12)*(s2-l24)*(s2-l14)); s3 = (l23 + l34 + l24) * 0.5; s3 = sqrt(s3*(s3-l23)*(s3-l34)*(s3-l24)); s4 = (l13 + l34 + l24) * 0.5; s4 = sqrt(s4*(s4-l13)*(s4-l34)*(s4-l14)); Real volume = Math::tetrahedron_volume(coord1,coord2,coord3,coord4); return 3*volume/(s1+s2+s3+s4); } /* -------------------------------------------------------------------------- */ inline void Math::barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter) { memset(barycenter, 0, spatial_dimension * sizeof(Real)); for (UInt n = 0; n < nb_points; ++n) { UInt offset = n * spatial_dimension; for (UInt i = 0; i < spatial_dimension; ++i) { barycenter[i] += coord[offset + i] / (Real) nb_points; } } } /* -------------------------------------------------------------------------- */ inline void Math::vector_2d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; } /* -------------------------------------------------------------------------- */ inline void Math::vector_3d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; - res[3] = y[2]-x[2]; + res[2] = y[2]-x[2]; } diff --git a/model/solid_mechanics/contact/contact_search_3d_explicit.cc b/model/solid_mechanics/contact/contact_search_3d_explicit.cc index e8db2fb06..8513edbed 100644 --- a/model/solid_mechanics/contact/contact_search_3d_explicit.cc +++ b/model/solid_mechanics/contact/contact_search_3d_explicit.cc @@ -1,468 +1,475 @@ /** * @file contact_search_3d_explicit.cc * @author David Kammer * @date Tue Oct 26 18:49:04 2010 * * @brief Specialization of the contact search structure for 3D within an explicit time scheme * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "regular_grid_neighbor_structure.hh" #include "contact_search_3d_explicit.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ ContactSearch3dExplicit::ContactSearch3dExplicit(Contact & contact, const ContactNeighborStructureType & neighbors_structure_type, const ContactSearchType & type, const ContactSearchID & id) : ContactSearch(contact, neighbors_structure_type, type, id), spatial_dimension(contact.getModel().getSpatialDimension()), mesh(contact.getModel().getFEM().getMesh()) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearch3dExplicit::findPenetration(const Surface & master_surface, PenetrationList & penetration_list) { AKANTU_DEBUG_IN(); /// get the NodesNeighborList for the given master surface std::map::iterator it_surface; for (it_surface = neighbors_structure.begin(); it_surface != neighbors_structure.end(); ++it_surface) { if(it_surface->first == master_surface) { break; } } AKANTU_DEBUG_ASSERT(it_surface != neighbors_structure.end(), "Master surface not found in this search object " << id); const NodesNeighborList & neighbor_list = dynamic_cast(it_surface->second->getNeighborList()); UInt nb_impactor_nodes = neighbor_list.impactor_nodes.getSize(); Vector * closest_master_nodes = new Vector(nb_impactor_nodes, 1); findClosestMasterNodes(master_surface, closest_master_nodes); UInt * closest_master_nodes_val = closest_master_nodes->values; /// get list of impactor and master nodes from neighbor list UInt * impactor_nodes_val = neighbor_list.impactor_nodes.values; //const Mesh & mesh = contact.getModel().getFEM().getMesh(); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; /// find existing surface element types UInt nb_types = type_list.size(); UInt nb_facet_types = 0; ElementType facet_type[_max_element_type]; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) == spatial_dimension) { ElementType current_facet_type = mesh.getFacetElementType(type); facet_type[nb_facet_types++] = current_facet_type; /// initialization of penetration list std::stringstream sstr_facets_offset; sstr_facets_offset << id << ":penetrated_facets_offset:" << current_facet_type; penetration_list.penetrated_facets_offset[current_facet_type] = new Vector(0, 1, sstr_facets_offset.str()); std::stringstream sstr_facets; sstr_facets << id << ":penetrated_facets:" << current_facet_type; penetration_list.penetrated_facets[current_facet_type] = new Vector(0, 1, sstr_facets.str()); std::stringstream sstr_normals; sstr_normals << id << ":facets_normals:" << current_facet_type; penetration_list.facets_normals[current_facet_type] = new Vector(0, spatial_dimension, sstr_normals.str()); std::stringstream sstr_gaps; sstr_gaps << id << ":gaps:" << current_facet_type; penetration_list.gaps[current_facet_type] = new Vector(0, 1, sstr_gaps.str()); std::stringstream sstr_projected_positions; sstr_projected_positions << id << ":projected_positions:" << current_facet_type; penetration_list.projected_positions[current_facet_type] = new Vector(0, spatial_dimension, sstr_projected_positions.str()); } } for(UInt in = 0; in < nb_impactor_nodes; ++in) { UInt current_impactor_node = impactor_nodes_val[in]; UInt closest_master_node = closest_master_nodes_val[in]; std::vector surface_elements; Vector * are_inside = new Vector(0, 1); Vector * are_in_projection_area = new Vector(0, 1); Element considered_element; for (UInt el_type = 0; el_type < nb_facet_types; ++el_type) { ElementType type = facet_type[el_type]; UInt * surface_id_val = mesh.getSurfaceId(type).values; const Vector & node_to_elements_offset = contact.getNodeToElementsOffset(type); const Vector & node_to_elements = contact.getNodeToElements(type); UInt * node_to_elements_offset_val = node_to_elements_offset.values; UInt * node_to_elements_val = node_to_elements.values; UInt min_element_offset = node_to_elements_offset_val[closest_master_node]; UInt max_element_offset = node_to_elements_offset_val[closest_master_node + 1]; considered_element.type = type; for(UInt el = min_element_offset; el < max_element_offset; ++el) { UInt surface_element = node_to_elements_val[el]; if(surface_id_val[surface_element] == master_surface) { bool is_inside; bool is_in_projection_area; checkPenetrationSituation(current_impactor_node, surface_element, type, is_inside, is_in_projection_area); considered_element.element = surface_element; surface_elements.push_back(considered_element); are_inside->push_back(is_inside); are_in_projection_area->push_back(is_in_projection_area); } } } UInt nb_penetrated_elements = 0; UInt nb_elements_type[_max_element_type]; memset(nb_elements_type, 0, sizeof(UInt) * _max_element_type); bool * are_inside_val = are_inside->values; bool * are_in_projection_area_val = are_in_projection_area->values; bool is_outside = false; UInt nb_surface_elements = static_cast(surface_elements.size()); // add all elements for which impactor is inside and in projection area for(UInt el = 0; el < nb_surface_elements; ++el) { if(are_inside_val[el] && are_in_projection_area_val[el]) { ElementType current_type = surface_elements.at(el).type; UInt current_element = surface_elements.at(el).element; penetration_list.penetrated_facets[current_type]->push_back(current_element); - Real normal[spatial_dimension]; - Real projected_position[spatial_dimension]; + Real normal[3]; + Real projected_position[3]; Real gap; computeComponentsOfProjection(current_impactor_node, current_element, current_type, normal, gap, projected_position); - penetration_list.facets_normals[current_element]->push_back(normal); - penetration_list.projected_positions[current_element]->push_back(projected_position); - penetration_list.gaps[current_element]->push_back(gap); + penetration_list.facets_normals[current_type]->push_back(normal); + penetration_list.projected_positions[current_type]->push_back(projected_position); + penetration_list.gaps[current_type]->push_back(gap); nb_penetrated_elements++; nb_elements_type[current_type]++; } } if(nb_penetrated_elements > 0) { for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) == spatial_dimension) { + penetration_list.penetrating_nodes.push_back(current_impactor_node); ElementType current_facet_type = mesh.getFacetElementType(type); penetration_list.penetrated_facets_offset[current_facet_type]->push_back(nb_elements_type[current_facet_type]); } } } // if there was no element which is inside and in projection area // check if node is not definitly outside else { for(UInt el = 0; el < nb_surface_elements && !is_outside; ++el) { if(!are_inside_val[el] && are_in_projection_area_val[el]) { is_outside = true; } - } - } - - // it is not definitly outside take all elements to which it is at least inside - if(!is_outside) { - for(UInt el = 0; el < nb_surface_elements; ++el) { - if(are_inside_val[el] && !are_in_projection_area_val[el]) { - - ElementType current_type = surface_elements.at(el).type; - UInt current_element = surface_elements.at(el).element; - penetration_list.penetrated_facets[current_type]->push_back(current_element); - - Real normal[spatial_dimension]; - Real projected_position[spatial_dimension]; - Real gap; - computeComponentsOfProjection(current_impactor_node, - current_element, - current_type, - normal, - gap, - projected_position); - - penetration_list.facets_normals[current_element]->push_back(normal); - penetration_list.projected_positions[current_element]->push_back(projected_position); - penetration_list.gaps[current_element]->push_back(gap); - - nb_penetrated_elements++; - nb_elements_type[current_type]++; + } + + // it is not definitly outside take all elements to which it is at least inside + if(!is_outside) { + bool found_inside_node = false; + for(UInt el = 0; el < nb_surface_elements; ++el) { + if(are_inside_val[el] && !are_in_projection_area_val[el]) { + + ElementType current_type = surface_elements.at(el).type; + UInt current_element = surface_elements.at(el).element; + penetration_list.penetrated_facets[current_type]->push_back(current_element); + + Real normal[3]; + Real projected_position[3]; + Real gap; + computeComponentsOfProjection(current_impactor_node, + current_element, + current_type, + normal, + gap, + projected_position); + + penetration_list.facets_normals[current_type]->push_back(normal); + penetration_list.projected_positions[current_type]->push_back(projected_position); + penetration_list.gaps[current_type]->push_back(gap); + + nb_penetrated_elements++; + nb_elements_type[current_type]++; + found_inside_node = true; + } } + if(found_inside_node) { + for(it = type_list.begin(); it != type_list.end(); ++it) { + ElementType type = *it; + if(mesh.getSpatialDimension(type) == spatial_dimension) { + penetration_list.penetrating_nodes.push_back(current_impactor_node); + ElementType current_facet_type = mesh.getFacetElementType(type); + penetration_list.penetrated_facets_offset[current_facet_type]->push_back(nb_elements_type[current_facet_type]); + } + } + } } - for(it = type_list.begin(); it != type_list.end(); ++it) { - ElementType type = *it; - if(mesh.getSpatialDimension(type) == spatial_dimension) { - ElementType current_facet_type = mesh.getFacetElementType(type); - penetration_list.penetrated_facets_offset[current_facet_type]->push_back(nb_elements_type[current_facet_type]); - } - } } delete are_in_projection_area; delete are_inside; } // make the offset table a real offset table for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) == spatial_dimension) { ElementType current_facet_type = mesh.getFacetElementType(type); UInt tmp_nb_facets = penetration_list.penetrated_facets_offset[current_facet_type]->getSize(); penetration_list.penetrated_facets_offset[current_facet_type]->resize(tmp_nb_facets+1); Vector & tmp_penetrated_facets_offset = *(penetration_list.penetrated_facets_offset[current_facet_type]); UInt * tmp_penetrated_facets_offset_val = tmp_penetrated_facets_offset.values; for (UInt i = 1; i < tmp_nb_facets; ++i) tmp_penetrated_facets_offset_val[i] += tmp_penetrated_facets_offset_val[i - 1]; for (UInt i = tmp_nb_facets; i > 0; --i) tmp_penetrated_facets_offset_val[i] = tmp_penetrated_facets_offset_val[i - 1]; tmp_penetrated_facets_offset_val[0] = 0; } } delete closest_master_nodes; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearch3dExplicit::findClosestMasterNodes(const Surface & master_surface, Vector * closest_master_nodes) { AKANTU_DEBUG_IN(); UInt * closest_master_nodes_val = closest_master_nodes->values; /// get the NodesNeighborList for the given master surface std::map::iterator it; for (it = neighbors_structure.begin(); it != neighbors_structure.end(); ++it) { if(it->first == master_surface) { break; } } AKANTU_DEBUG_ASSERT(it != neighbors_structure.end(), "Master surface not found in this search object " << id); const NodesNeighborList & neighbor_list = dynamic_cast(it->second->getNeighborList()); /// get list of impactor and master nodes from neighbor list UInt * impactor_nodes_val = neighbor_list.impactor_nodes.values; UInt * master_nodes_offset_val = neighbor_list.master_nodes_offset.values; UInt * master_nodes_val = neighbor_list.master_nodes.values; /// loop over all impactor nodes and find for each the closest master node UInt nb_impactor_nodes = neighbor_list.impactor_nodes.getSize(); for(UInt imp = 0; imp < nb_impactor_nodes; ++imp) { UInt current_impactor_node = impactor_nodes_val[imp]; UInt min_offset = master_nodes_offset_val[imp]; UInt max_offset = master_nodes_offset_val[imp + 1]; Real min_square_distance = std::numeric_limits::max(); UInt closest_master_node; for(UInt mn = min_offset; mn < max_offset; ++mn) { UInt current_master_node = master_nodes_val[mn]; Real square_distance = computeSquareDistanceBetweenNodes(current_impactor_node, current_master_node); if(min_square_distance > square_distance) { min_square_distance = square_distance; closest_master_node = current_master_node; } } closest_master_nodes_val[imp] = closest_master_node; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearch3dExplicit::computeComponentsOfProjection(const UInt impactor_node, const UInt surface_element, const ElementType type, Real * normal, Real & gap, Real * projected_position) { AKANTU_DEBUG_IN(); switch(type) { case _triangle_3: { computeComponentsOfProjectionTriangle3(impactor_node, surface_element, normal, gap, projected_position); break; } case _not_defined: { AKANTU_DEBUG_ERROR("Not a valid surface element type : " << type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearch3dExplicit::checkPenetrationSituation(const UInt impactor_node, const UInt surface_element, const ElementType type, bool & is_inside, bool & is_in_projection_area) { AKANTU_DEBUG_IN(); switch(type) { case _triangle_3: { checkPenetrationSituationTriangle3(impactor_node, surface_element, is_inside, is_in_projection_area); break; } case _not_defined: { AKANTU_DEBUG_ERROR("Not a valid surface element type : " << type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearch3dExplicit::computeComponentsOfProjectionTriangle3(const UInt impactor_node, const UInt surface_element, Real * normal, Real & gap, Real * projected_position) { AKANTU_DEBUG_IN(); const UInt dim = spatial_dimension; const ElementType type = _triangle_3; const UInt nb_nodes_element = Mesh::getNbNodesPerElement(type); Real * current_position = contact.getModel().getCurrentPosition().values; UInt * connectivity = mesh.getConnectivity(type).values; UInt node_1 = surface_element * nb_nodes_element; Real * position_node_1 = &(current_position[connectivity[node_1 + 0] * spatial_dimension]); Real * position_node_2 = &(current_position[connectivity[node_1 + 1] * spatial_dimension]); Real * position_node_3 = &(current_position[connectivity[node_1 + 2] * spatial_dimension]); Real * position_impactor_node = &(current_position[impactor_node * spatial_dimension]); /// compute the normal of the face - Real vector_1[dim]; - Real vector_2[dim]; - + Real vector_1[3]; + Real vector_2[3]; Math::vector_3d(position_node_1, position_node_2, vector_1); /// @todo: check if correct order of nodes !! Math::vector_3d(position_node_1, position_node_3, vector_2); Math::normal3(vector_1, vector_2, normal); /// compute the gap between impactor and face /// gap positive if impactor outside of surface - Real vector_node_1_impactor[dim]; + Real vector_node_1_impactor[3]; Math::vector_3d(position_node_1, position_impactor_node, vector_node_1_impactor); gap = Math::vectorDot3(vector_node_1_impactor, normal); /// compute the projected position of the impactor node onto the face for(UInt i=0; i < dim; ++i) { projected_position[i] = position_impactor_node[i] - gap * normal[i]; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactSearch3dExplicit::checkPenetrationSituationTriangle3(const UInt impactor_node, const UInt surface_element, bool & is_inside, bool & is_in_projection_area) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(spatial_dimension == 3, "wrong spatial dimension (=" << spatial_dimension << ") for checkPenetrationSituationTriangle3"); const UInt dim = spatial_dimension; const ElementType type = _triangle_3; const UInt nb_nodes_element = Mesh::getNbNodesPerElement(type); Real * current_position = contact.getModel().getCurrentPosition().values; UInt * connectivity = mesh.getConnectivity(type).values; - - Real normal[dim]; - Real gap; - Real projected_position[dim]; + + Real gap; + Real normal[3]; + Real projected_position[3]; computeComponentsOfProjectionTriangle3(impactor_node, surface_element, normal, gap, projected_position); // ------------------------------------------------------- /// Find if impactor node is inside or outside of the face // ------------------------------------------------------- if(gap < 0.0) is_inside = true; else is_inside = false; // ---------------------------------------------------- /// Find if impactor node is in projection area of face // ---------------------------------------------------- UInt node_1 = surface_element * nb_nodes_element; Real * position_node_1 = &(current_position[connectivity[node_1 + 0] * spatial_dimension]); Real * position_node_2 = &(current_position[connectivity[node_1 + 1] * spatial_dimension]); Real * position_node_3 = &(current_position[connectivity[node_1 + 2] * spatial_dimension]); Real triangle_area; Real test_area = 0.0; - Real tmp_vector_1[dim]; - Real tmp_vector_2[dim]; - Real tmp_vector_3[dim]; + Real tmp_vector_1[3]; + Real tmp_vector_2[3]; + Real tmp_vector_3[3]; // find area of triangle Math::vector_3d(position_node_1, position_node_2, tmp_vector_1); Math::vector_3d(position_node_1, position_node_3, tmp_vector_2); Math::vectorProduct3(tmp_vector_1, tmp_vector_2, tmp_vector_3); triangle_area = 0.5 * Math::norm3(tmp_vector_3); // compute areas of projected position and two nodes of master triangle + UInt nb_sub_areas = nb_nodes_element; UInt node_order[4]; node_order[0] = 0; node_order[1] = 1; node_order[2] = 2; node_order[3] = 0; - for(UInt i=0; i < spatial_dimension; ++i) { + for(UInt i=0; i < nb_sub_areas; ++i) { position_node_1 = &(current_position[connectivity[node_1 + node_order[i ]] * spatial_dimension]); position_node_2 = &(current_position[connectivity[node_1 + node_order[i+1]] * spatial_dimension]); Math::vector_3d(projected_position, position_node_1, tmp_vector_1); Math::vector_3d(projected_position, position_node_2, tmp_vector_2); Math::vectorProduct3(tmp_vector_1, tmp_vector_2, tmp_vector_3); test_area += 0.5 * Math::norm3(tmp_vector_3); } // the projection is outside if the test area is larger than the area of the triangle - if(test_area > triangle_area) + Real tolerance = std::numeric_limits::epsilon(); + if((test_area - triangle_area) > tolerance) is_in_projection_area = false; else is_in_projection_area = true; AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 42b9f1a49..560e2c9fe 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,41 +1,41 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # # # @section DESCRIPTION # #=============================================================================== INCLUDE(${AKANTU_CMAKE_DIR}/AkantuTestAndExamples.cmake) #=============================================================================== # List of tests #=============================================================================== add_akantu_test(test_vector "Test akantu vector") add_akantu_test(test_static_memory "Test static memory") add_akantu_test(test_fem "Test finite element functionalties") add_akantu_test(test_mesh_io "Test mesh io object") add_akantu_test(test_facet_extraction "Test mesh utils facet extraction") add_akantu_test(test_model "Test model objects") #add_akantu_test(test_solid_mechanics_model "Test solid mechanics object") #add_akantu_test(test_solver "Test solver function") if(AKANTU_BUILD_CONTACT) - add_akantu_test(test_contact_neighbor_structure "Test contact neighbor structure") + #add_akantu_test(test_contact_neighbor_structure "Test contact neighbor structure") add_akantu_test(test_surface_extraction "Test mesh utils surface extraction") endif(AKANTU_BUILD_CONTACT) if(AKANTU_SCOTCH_ON) add_akantu_test(test_mesh_partitionate "Test mesh partition creation") endif(AKANTU_SCOTCH_ON) if(AKANTU_MPI_ON) add_akantu_test(test_synchronizer "Test synchronizers") endif(AKANTU_MPI_ON) diff --git a/test/test_model/test_solid_mechanics_model/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/CMakeLists.txt index c5c721ce0..ae2a6f778 100644 --- a/test/test_model/test_solid_mechanics_model/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/CMakeLists.txt @@ -1,68 +1,72 @@ #=============================================================================== # @file CMakeLists.txt # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # # # @section DESCRIPTION # #=============================================================================== add_subdirectory("test_materials") +if(AKANTU_BUILD_CONTACT) + add_subdirectory("test_contact") +endif(AKANTU_BUILD_CONTACT) + #=============================================================================== add_mesh(test_solid_mechanics_model_square_mesh triangle.geo 2 1) add_mesh(test_solid_mechanics_model_circle_mesh circle.geo 2 2) register_test(test_solid_mechanics_model_square test_solid_mechanics_model_square.cc) add_dependencies(test_solid_mechanics_model_square test_solid_mechanics_model_square_mesh test_solid_mechanics_model_circle_mesh) register_test(test_solid_mechanics_model_circle_2 test_solid_mechanics_model_circle_2.cc) #=============================================================================== add_mesh(test_bar_traction_2d_mesh1 bar.geo 2 1 OUTPUT bar1.msh) add_mesh(test_bar_traction_2d_mesh2 bar.geo 2 2 OUTPUT bar2.msh) register_test(test_solid_mechanics_model_bar_traction2d test_solid_mechanics_model_bar_traction2d.cc) add_dependencies(test_solid_mechanics_model_bar_traction2d test_bar_traction_2d_mesh1 test_bar_traction_2d_mesh2) if(AKANTU_SCOTCH_ON AND AKANTU_MPI_ON) register_test(test_solid_mechanics_model_bar_traction2d_parallel test_solid_mechanics_model_bar_traction2d_parallel.cc) add_dependencies(test_solid_mechanics_model_bar_traction2d_parallel test_bar_traction_2d_mesh2) add_mesh(test_solid_mechanics_model_segment_mesh segment.geo 1 2) register_test(test_solid_mechanics_model_segment_parallel test_solid_mechanics_model_segment_parallel.cc) add_dependencies(test_solid_mechanics_model_segment_parallel test_solid_mechanics_model_segment_mesh) endif() #=============================================================================== add_mesh(test_cube3d_mesh1 cube.geo 3 1 OUTPUT cube1.msh) add_mesh(test_cube3d_mesh2 cube.geo 3 2 OUTPUT cube2.msh) register_test(test_solid_mechanics_model_cube3d test_solid_mechanics_model_cube3d.cc) add_dependencies(test_solid_mechanics_model_cube3d test_cube3d_mesh1) register_test(test_solid_mechanics_model_cube3d_tetra10 test_solid_mechanics_model_cube3d_tetra10.cc) add_dependencies(test_solid_mechanics_model_cube3d_tetra10 test_cube3d_mesh2) #=============================================================================== configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/material.dat ${CMAKE_CURRENT_BINARY_DIR}/material.dat COPYONLY ) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_model/test_solid_mechanics_model/test_contact/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_contact/CMakeLists.txt new file mode 100644 index 000000000..09b4b9777 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_contact/CMakeLists.txt @@ -0,0 +1,17 @@ +#=============================================================================== +# @file CMakeLists.txt +# @author David Kammer +# @date Fri Dec 03 11:09:59 2010 +# +# @section LICENSE +# +# +# +# @section DESCRIPTION +# +#=============================================================================== + +#=============================================================================== + +add_subdirectory("test_neighbor_structure") +add_subdirectory("test_contact_search") diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/CMakeLists.txt new file mode 100644 index 000000000..958db93f8 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/CMakeLists.txt @@ -0,0 +1,39 @@ +#=============================================================================== +# @file CMakeLists.txt +# @author David Kammer +# @date Fri Dec 03 12:02:24 2010 +# +# @section LICENSE +# +# +# +# @section DESCRIPTION +# +#=============================================================================== + +#=============================================================================== +#add_mesh(test_regular_grid_2d_mesh squares.geo 2 1) + +#register_test(test_regular_grid_triangle_3 test_regular_grid_triangle_3.cc) +#add_dependencies(test_regular_grid_triangle_3 test_regular_grid_2d_mesh) + +#=============================================================================== +add_mesh(test_3d_explicit_mesh cubes.geo 3 1) + +register_test(test_3d_explicit_tetrahedron_4 test_3d_explicit_tetrahedron_4.cc) +add_dependencies(test_3d_explicit_tetrahedron_4 test_3d_explicit_mesh) + +#=============================================================================== +#add_mesh(test_regular_grid_3d_mesh cubes.geo 3 1) + +#register_test(test_regular_grid_tetrahedron_4_nodes test_regular_grid_tetrahedron_4_nodes.cc) +#add_dependencies(test_regular_grid_tetrahedron_4_nodes test_regular_grid_3d_mesh) + +#=============================================================================== +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/material.dat + ${CMAKE_CURRENT_BINARY_DIR}/material.dat + COPYONLY + ) + +file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_contact_neighbor_structure/cubes.geo.small b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/cubes.geo similarity index 74% copy from test/test_contact_neighbor_structure/cubes.geo.small copy to test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/cubes.geo index b65d24a12..2939e7cc2 100644 --- a/test/test_contact_neighbor_structure/cubes.geo.small +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/cubes.geo @@ -1,74 +1,81 @@ -h = 1; +h = 0.05; +h1 = 0.4; -Point(1) = {0, 0, 0, h}; -Point(2) = {1, 0, 0, h}; -Point(3) = {1, 1, 0, h}; -Point(4) = {0, 1, 0, h}; +Point(1) = {0, 0, 0.7, h}; +Point(2) = {1, 0, 0.7, h}; +Point(3) = {1, 1, 0.7, h}; +Point(4) = {0, 1, 0.7, h}; Point(5) = {0, 0, 1, h}; Point(6) = {1, 0, 1, h}; Point(7) = {1, 1, 1, h}; Point(8) = {0, 1, 1, h}; -Point(9) = {1.01+0, 0, 0, h}; -Point(10) = {1.01+1, 0, 0, h}; -Point(11) = {1.01+1, 1, 0, h}; -Point(12) = {1.01+0, 1, 0, h}; -Point(13) = {1.01+0, 0, 1, h}; -Point(14) = {1.01+1, 0, 1, h}; -Point(15) = {1.01+1, 1, 1, h}; -Point(16) = {1.01+0, 1, 1, h}; +Point(9) = {0.45, 0.45, 1.001, h1}; +Point(10) = {0.55, 0.45, 1.001, h1}; +Point(11) = {0.55, 0.55, 1.001, h1}; +Point(12) = {0.45, 0.55, 1.001, h1}; +Point(13) = {0.45, 0.45, 1.101, h1}; +Point(14) = {0.55, 0.45, 1.101, h1}; +Point(15) = {0.55, 0.55, 1.101, h1}; +Point(16) = {0.45, 0.55, 1.101, h1}; Line(1) = {5, 6}; Line(2) = {6, 7}; Line(3) = {7, 8}; Line(4) = {8, 5}; Line(5) = {1, 2}; Line(6) = {2, 3}; Line(7) = {3, 4}; Line(8) = {4, 1}; Line(9) = {5, 1}; Line(10) = {6, 2}; Line(11) = {7, 3}; Line(12) = {8, 4}; + Line(13) = {13, 14}; Line(14) = {14, 15}; Line(15) = {15, 16}; Line(16) = {16, 13}; Line(17) = {9, 10}; Line(18) = {10, 11}; Line(19) = {11, 12}; Line(20) = {12, 9}; Line(21) = {13, 9}; Line(22) = {14, 10}; Line(23) = {15, 11}; Line(24) = {16, 12}; - Line Loop(25) = {4, 1, 2, 3}; Plane Surface(26) = {25}; Line Loop(27) = {7, 8, 5, 6}; Plane Surface(28) = {27}; Line Loop(29) = {12, -7, -11, 3}; Plane Surface(30) = {29}; Line Loop(31) = {12, 8, -9, -4}; Plane Surface(32) = {31}; Line Loop(33) = {1, 10, -5, -9}; Plane Surface(34) = {33}; Line Loop(35) = {2, 11, -6, -10}; Plane Surface(36) = {35}; + Line Loop(37) = {16, 13, 14, 15}; Plane Surface(38) = {37}; Line Loop(39) = {20, 17, 18, 19}; Plane Surface(40) = {39}; Line Loop(41) = {24, 20, -21, -16}; Plane Surface(42) = {41}; Line Loop(43) = {13, 22, -17, -21}; Plane Surface(44) = {43}; Line Loop(45) = {14, 23, -18, -22}; Plane Surface(46) = {45}; Line Loop(47) = {15, 24, -19, -23}; Plane Surface(48) = {47}; + Surface Loop(49) = {32, 30, 28, 34, 26, 36}; Volume(50) = {49}; + Surface Loop(51) = {42, 48, 38, 44, 46, 40}; Volume(52) = {51}; + +Physical Volume(0) = {50}; +Physical Volume(1) = {52}; diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/material.dat b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/material.dat new file mode 100644 index 000000000..0876b3860 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/material.dat @@ -0,0 +1,13 @@ +material elastic [ + name = steel_1 + rho = 7800 # density + E = 2.1e11 # young's modulus + nu = 0.3 # poisson's ratio +] + +material elastic [ + name = steel_2 + rho = 7800 # density + E = 2.1e11 # young's modulus + nu = 0.3 # poisson's ratio +] diff --git a/test/test_contact_neighbor_structure/test_contact_regular_grid_tetrahedron_4_nodes.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_3d_explicit_tetrahedron_4.cc similarity index 52% copy from test/test_contact_neighbor_structure/test_contact_regular_grid_tetrahedron_4_nodes.cc copy to test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_3d_explicit_tetrahedron_4.cc index 47977239f..3d0d5a97e 100644 --- a/test/test_contact_neighbor_structure/test_contact_regular_grid_tetrahedron_4_nodes.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_3d_explicit_tetrahedron_4.cc @@ -1,126 +1,189 @@ /** * @file test_contact_regular_grid.cc * @author David Kammer - * @date Tue Oct 26 16:58:42 2010 + * @date Fri Dec 03 12:11:42 2010 * - * @brief test regular grid neighbor structure for 3d case + * @brief test contact search for 3d case in explicit * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" - - +#include "contact_search.hh" +#include "contact_search_3d_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; + const ElementType element_type = _tetrahedron_4; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("cubes.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); - + + UInt max_steps = 2; unsigned int nb_nodes = my_mesh.getNbNodes(); - /// dump facet and surface information to paraview #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "tetrahedron_4_nodes_test-surface-extraction"); dumper.SetConnectivity((int*)my_mesh.getConnectivity(_tetrahedron_4).values, TETRA1, my_mesh.getNbElement(_tetrahedron_4), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors - memset(my_model.getForce().values, 0, 3*nb_nodes*sizeof(Real)); - memset(my_model.getVelocity().values, 0, 3*nb_nodes*sizeof(Real)); - memset(my_model.getAcceleration().values, 0, 3*nb_nodes*sizeof(Real)); - memset(my_model.getDisplacement().values, 0, 3*nb_nodes*sizeof(Real)); + memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); + memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); + memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); + memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); + + Real * displacement = my_model.getDisplacement().values; my_model.readMaterials("material.dat"); my_model.initMaterials(); my_model.initModel(); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_3d_expli, _cst_3d_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 0; - + Surface impactor = 1; my_contact->addMasterSurface(master); - my_contact->initNeighborStructure(master); - + const NodesNeighborList & my_neighbor_list = dynamic_cast(my_contact->getContactSearch().getContactNeighborStructure(master).getNeighborList()); UInt nb_nodes_neigh = my_neighbor_list.impactor_nodes.getSize(); Vector impact_nodes = my_neighbor_list.impactor_nodes; UInt * impact_nodes_val = impact_nodes.values; /// print impactor nodes std::cout << "we have " << nb_nodes_neigh << " impactor nodes:"; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " " << impact_nodes_val[i]; } std::cout << std::endl; UInt * master_nodes_offset_val = my_neighbor_list.master_nodes_offset.values; UInt * master_nodes_val = my_neighbor_list.master_nodes.values; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " Impactor node: " << impact_nodes_val[i] << " has master nodes:"; for(UInt mn = master_nodes_offset_val[i]; mn < master_nodes_offset_val[i+1]; ++mn) { std::cout << " " << master_nodes_val[mn]; } std::cout << std::endl; - } + } + + my_contact->initSearch(); // does nothing so far + + std::cout << std::endl << "epsilon = " << std::numeric_limits::epsilon() << std::endl; + + /* ------------------------------------------------------------------------ */ + /* Main loop */ + /* ------------------------------------------------------------------------ */ + for(UInt s = 1; s <= max_steps; ++s) { + + std::cout << std::endl << "passing step " << s << "/" << max_steps << std::endl; + + /// apply a displacement to the slave body + if(s == 2) { + Real * coord = my_mesh.getNodes().values; + for(UInt n = 0; n < nb_nodes; ++n) { + if(coord[n*dim + 2] > 1.0) { + displacement[n*dim+2] = -0.01; + } + } + /* + UInt nb_elements = my_mesh.getNbElement(element_type); + UInt nb_nodes_element = my_mesh.getNbNodesPerElement(element_type); + Vector element_mat = my_model.getElementMaterial(element_type); + UInt * element_mat_val = element_mat.values; + UInt * connectivity = my_mesh.getConnectivity(element_type).values; + for(UInt el = 0; el < nb_elements; ++el) { + std::cout << "element: " << el << " with mat: " << element_mat_val[el] << std::endl; + if(element_mat_val[el] == impactor) { + for(UInt n = 0; n < nb_nodes_element; ++n) { + displacement[connectivity[el * nb_nodes_element + n]+2] = -0.2; + } + } + }*/ + } + + /// central difference predictor + my_model.explicitPred(); + + /// compute the residual + my_model.updateResidual(); + + /// compute the penetration list + PenetrationList * my_penetration_list = new PenetrationList(); + const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); + + UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); + Vector pen_nodes = my_penetration_list->penetrating_nodes; + UInt * pen_nodes_val = pen_nodes.values; + std::cout << "we have " << nb_nodes_pen << " penetrating nodes:"; + for (UInt i = 0; i < nb_nodes_pen; ++i) + std::cout << " " << pen_nodes_val[i]; + std::cout << std::endl; + delete my_penetration_list; + + /// compute the acceleration + my_model.updateAcceleration(); + + /// central difference corrector + my_model.explicitCorr(); + } delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_contact_neighbor_structure/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/CMakeLists.txt similarity index 50% rename from test/test_contact_neighbor_structure/CMakeLists.txt rename to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/CMakeLists.txt index 1b7a31bad..756082bc0 100644 --- a/test/test_contact_neighbor_structure/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/CMakeLists.txt @@ -1,30 +1,39 @@ #=============================================================================== # @file CMakeLists.txt # @author David Kammer # @date Tue Oct 26 16:40:24 2010 # # @section LICENSE # # # # @section DESCRIPTION # #=============================================================================== #=============================================================================== add_mesh(test_regular_grid_2d_mesh squares.geo 2 1) -register_test(test_contact_regular_grid_triangle_3 test_contact_regular_grid_triangle_3.cc) -add_dependencies(test_contact_regular_grid_triangle_3 test_regular_grid_2d_mesh) +register_test(test_regular_grid_triangle_3 test_regular_grid_triangle_3.cc) +add_dependencies(test_regular_grid_triangle_3 test_regular_grid_2d_mesh) #=============================================================================== add_mesh(test_regular_grid_3d_mesh cubes.geo 3 1) -register_test(test_contact_regular_grid_tetrahedron_4 test_contact_regular_grid_tetrahedron_4.cc) -add_dependencies(test_contact_regular_grid_tetrahedron_4 test_regular_grid_3d_mesh) +register_test(test_regular_grid_tetrahedron_4 test_regular_grid_tetrahedron_4.cc) +add_dependencies(test_regular_grid_tetrahedron_4 test_regular_grid_3d_mesh) #=============================================================================== #add_mesh(test_regular_grid_3d_mesh cubes.geo 3 1) -register_test(test_contact_regular_grid_tetrahedron_4_nodes test_contact_regular_grid_tetrahedron_4_nodes.cc) -add_dependencies(test_contact_regular_grid_tetrahedron_4_nodes test_regular_grid_3d_mesh) \ No newline at end of file +register_test(test_regular_grid_tetrahedron_4_nodes test_regular_grid_tetrahedron_4_nodes.cc) +add_dependencies(test_regular_grid_tetrahedron_4_nodes test_regular_grid_3d_mesh) + +#=============================================================================== +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/material.dat + ${CMAKE_CURRENT_BINARY_DIR}/material.dat + COPYONLY + ) + +file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_contact_neighbor_structure/cubes.geo b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/cubes.geo similarity index 100% rename from test/test_contact_neighbor_structure/cubes.geo rename to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/cubes.geo diff --git a/test/test_contact_neighbor_structure/cubes.geo.large b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/cubes.geo.large similarity index 100% rename from test/test_contact_neighbor_structure/cubes.geo.large rename to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/cubes.geo.large diff --git a/test/test_contact_neighbor_structure/cubes.geo.small b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/cubes.geo.small similarity index 100% rename from test/test_contact_neighbor_structure/cubes.geo.small rename to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/cubes.geo.small diff --git a/test/test_contact_neighbor_structure/material.dat b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/material.dat similarity index 100% rename from test/test_contact_neighbor_structure/material.dat rename to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/material.dat diff --git a/test/test_contact_neighbor_structure/squares.geo b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/squares.geo similarity index 100% rename from test/test_contact_neighbor_structure/squares.geo rename to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/squares.geo diff --git a/test/test_contact_neighbor_structure/test_contact_regular_grid_tetrahedron_4.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_tetrahedron_4.cc similarity index 100% rename from test/test_contact_neighbor_structure/test_contact_regular_grid_tetrahedron_4.cc rename to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_tetrahedron_4.cc diff --git a/test/test_contact_neighbor_structure/test_contact_regular_grid_tetrahedron_4_nodes.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_tetrahedron_4_nodes.cc similarity index 100% rename from test/test_contact_neighbor_structure/test_contact_regular_grid_tetrahedron_4_nodes.cc rename to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_tetrahedron_4_nodes.cc diff --git a/test/test_contact_neighbor_structure/test_contact_regular_grid_triangle_3.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_triangle_3.cc similarity index 100% rename from test/test_contact_neighbor_structure/test_contact_regular_grid_triangle_3.cc rename to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_triangle_3.cc