contact_searching.cpp

vs.lib is a math library in C++ with a set of linear algebra and integrable / differentiable objects

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <iomanip.h>
#include "vs_h.h"
#include "dynamic_array.h"
#include "omega_h.h"
#include "u_h.h"
#include "element_formulation.h"
#include "matrix_representation.h"
#include "omega_h_n.h"
#include "global_discretization_gamma_h_n.h"
#include "contact_searching.h"

Contact_Omega_eh::Contact_Omega_eh(int* ae, double* pp, double* nv, Omega_eh& oh)
	: Omega_eh(oh) {
   	for(int i = 0; i < 2; i++) the_associate_element[i] = ae[i];
      for(i = 0; i < 4; i++) the_projection_parameter[i] = pp[i];
      for(i = 0; i < 4; i++) the_normal_vector[i] = nv[i];
}

int setup_order_table(int *order_table, int size,
	Global_Discretization_Gamma_h_i &a_gd, Global_Discretization_Gamma_h_i &b_gd, double tolerance_gap = 1.e-3);

void define_contact_node(Contact_Omega_h_i &interface_oh, int *order_table, int size, int order_table_count,
	Global_Discretization_Gamma_h_i &a_gd, Global_Discretization_Gamma_h_i &b_gd, double tolerance_gap = 1.e-3);

void define_contact_element(Contact_Omega_h_i &interface_oh, int *order_table, int size, int order_table_count,
	Global_Discretization_Gamma_h_i &a_gd, Global_Discretization_Gamma_h_i &b_gd, double tolerance_gap = 1.e-3);
   
void Contact_Omega_h_i::contact_searching_algorithm(Contact_Omega_h_i &interface_oh, Global_Discretization_Gamma_h_i &a_gd,
	Global_Discretization_Gamma_h_i &b_gd, U_h& lambda_h, gh_on_Gamma_h_i& interface_gh, double tolerance_gap) {
	// contact searching
   Gamma_h_i &a = a_gd.gamma_h_i(), &b = b_gd.gamma_h_i();
	int a_total_node_no = a.total_element_no()+1,
       b_total_node_no = b.total_element_no()+1,
       max_size = a_total_node_no+b_total_node_no;
   // instantiation; allocate memory
   int **order_table = new int*[2];
   order_table[0] = new int[max_size*2];
   order_table[1] = order_table[0] + max_size;
   for(int i = 0; i < 2; i++) // initialization; define default value
   for(int j = 0; j < max_size; j++) order_table[i][j] = -1;
   int order_table_count = setup_order_table(order_table[0], max_size, a_gd, b_gd, tolerance_gap);

   // define contact nodes
   define_contact_node(interface_oh, order_table[0], max_size, order_table_count, a_gd, b_gd, tolerance_gap);

   // define contact elements
   define_contact_element(interface_oh, order_table[0], max_size, order_table_count, a_gd, b_gd, tolerance_gap);
   
	extern ofstream ofs;
	ofs << "order_table size: " << max_size << endl;
	for(i = 0; i < 2; i++) {
		for(int j = 0; j < max_size; j++) ofs << order_table[i][j] << ", ";
		ofs << endl;
	}
   delete [] order_table[0];
   delete [] order_table;

   // re-define lambda_h; the lagrangian variables, and its B.C.
   lambda_h.u_h_array().~Dynamic_Array<Nodal_Value>(); // explicit call to destructor
   lambda_h.__initialization(2, interface_oh);
   interface_gh.gh_array().~Dynamic_Array<Nodal_Constraint>();
   interface_gh.__initialization(2, interface_oh);
}
   
void Contact_Omega_h_i::contact_searching_algorithm(Contact_Omega_h_i &interface_oh, Global_Discretization_Gamma_h_i &a_gd,
	Global_Discretization_Gamma_h_i &b_gd, U_h& lambda_h, gh_on_Gamma_h_i& interface_gh,
   Matrix_Representation& mr_K_0, Matrix_Representation& mr_K_1,
   Matrix_Representation_Couple& mrc_Q_0, Matrix_Representation_Couple& mrc_Q_1,
   double tolerance_gap) {

   contact_searching_algorithm(interface_oh, a_gd, b_gd, lambda_h, interface_gh);
   
   // rhs of mr needs to be cleaned
   mr_K_0.~Matrix_Representation();
   mr_K_0.__initialization(0);
   mr_K_1.~Matrix_Representation();
   mr_K_1.__initialization(0);

   
   mrc_Q_0.~Matrix_Representation_Couple();
   ((Matrix_Representation*)&mrc_Q_0)->__initialization("NONE_LHS_RHS");
   mrc_Q_0.__initialization(0, 0, &(mr_K_0.rhs()) );

   
   mrc_Q_1.~Matrix_Representation_Couple();
   ((Matrix_Representation*)&mrc_Q_1)->__initialization("NONE_LHS_RHS");
   mrc_Q_1.__initialization(0, &(mrc_Q_0.rhs()), &(mr_K_1.rhs()) );
}

//------------------------------------------------------------------------------
//
// Geometrical Utitilies for Contact Searching Algorithm
//
//------------------------------------------------------------------------------
// function prototypes of geometry utilities
int node_on_node(Node&, Global_Discretization_Gamma_h_i&, Node&,
	Global_Discretization_Gamma_h_i&, double = 1.e-3, C0& = C0(2,(double*)0));
int node_on_segment(Node&, Global_Discretization_Gamma_h_i&, Node&, Node&,
	Global_Discretization_Gamma_h_i&, double &alpha, double tolerance_gap = 1.e-3,
   C0 &ret_coord = C0(2,(double*)0));
int node_projection(Node &projectile_node, Global_Discretization_Gamma_h_i& projectile,
	Global_Discretization_Gamma_h_i &target_gd, double &alpha,
	double tolerance_gap = 1.e-3, int target_first_node_no = 0,
   C0 &ret_coord = C0(2,(double*)0));

enum Contact_Condition_Code {NODE_ON_NODE, A_NODE_LEADING, B_NODE_LEADING, PROJECTION_FAIL}; // contact condition codes
Contact_Condition_Code contact_node_searching(Global_Discretization_Gamma_h_i &a_gd, int &a_node_no,
	Global_Discretization_Gamma_h_i &b_gd, int &b_node_no, double tolerance_gap = 1.e-3) {
   Gamma_h_i &a = a_gd.gamma_h_i(), &b = b_gd.gamma_h_i();
   int a_total_node_no = a.total_element_no()+1,
       b_total_node_no = b.total_element_no()+1,
   	 a_segment, b_segment;
   Node a_node, b_node;
   do {
      if(a_node_no < a_total_node_no-1) a_node = a[a(a_node_no)[0]]; // not the last node
   	else a_node = a[a(a_node_no-1)[1]];
      double alpha;
   	if( (b_segment = node_projection(a_node, a_gd, b_gd, alpha, tolerance_gap, 0)) > -1) break;
   } while(b_segment == -1 && a_node_no++ < a_total_node_no); // skip to next node
   do {
      if(b_node_no < b_total_node_no-1) b_node = b[b(b_node_no)[0]]; // not the last node
      else b_node = b[b(b_node_no-1)[1]];
      double alpha;
   	if( (a_segment = node_projection(b_node, b_gd, a_gd, alpha, tolerance_gap, 0)) > -1) break;
   } while(a_segment == -1 && b_node_no++ < b_total_node_no); // skip to next node
   if(a_segment == -1 && b_segment == -1) return PROJECTION_FAIL;
   if(node_on_node(a_node, a_gd, b_node, b_gd, tolerance_gap, C0(2,(double*)0))) return NODE_ON_NODE;
   if((a_node_no <= a_segment) || a_segment == -1) return A_NODE_LEADING;
   return B_NODE_LEADING;
}

int node_projection(Node &projectile_node, Global_Discretization_Gamma_h_i& projectile_gd,
	Global_Discretization_Gamma_h_i &target_gd, double &alpha, double tolerance_gap,
	int target_first_node_no, C0 &ret_coord) {
   Gamma_h_i &target = target_gd.gamma_h_i();
   int target_count = target_first_node_no;
	do{
		Omega_eh &target_elem = target(target_count);
      Node &node1 = target[target_elem[0]];
      Node &node2 = target[target_elem[1]];
      if(node_on_segment(projectile_node, projectile_gd, node1, node2, target_gd,
      	alpha, tolerance_gap, ret_coord))  break;            // node on segment; 0 < alpha < 1
      else if(node_on_node(projectile_node, projectile_gd, node1, target_gd, tolerance_gap,
                           ret_coord)) { alpha = 0.0; break; } // node on first node
      else if(node_on_node(projectile_node, projectile_gd, node2, target_gd, tolerance_gap,
                           ret_coord)) { alpha = 1.0; break; } // node on second node
   } while(target_count++ < target.total_element_no());
   if((target_count == target.total_element_no()+1) &&
      (alpha < 0.0 || alpha > 1.0) )
      	return -1;                // projection failed
   return target_count;
}

int node_on_node(Node& a, Global_Discretization_Gamma_h_i& a_gd, Node& b,
	Global_Discretization_Gamma_h_i& b_gd, double tolerance_gap, C0 &ret_coord) {
   int a_node_no = a.node_no(),
       b_node_no = b.node_no();
	C0 x(2, a.value()),                            // initial coordinates
   	y(2, b.value()),
      xu(2, a_gd.u_h()[a_node_no].value()),       // displacements
      yu(2, b_gd.u_h()[b_node_no].value()),
      x_coord = x+xu,                             // lagrangian(material) coordinates
      y_coord = y+yu;
   if((double)norm(x_coord-y_coord, 2) <= tolerance_gap) {  // if close enough
   	ret_coord = (x_coord+y_coord)/2.0;          // return contact node coordinates
   	return TRUE;
   }
   return FALSE;
}

int node_on_segment(Node &c, Global_Discretization_Gamma_h_i& node_gd, Node &a, Node &b,
	Global_Discretization_Gamma_h_i& seg_gd, double &alpha, double tolerance_gap, C0 &ret_coord) {
	// examine if c can be projected onto segment a-b
   int a_node_no = a.node_no(),
       b_node_no = b.node_no(),
       c_node_no = c.node_no();
	C0 x(2, c.value()),                                  // initial coordinates
      x_1(2, a.value()),
      x_2(2, b.value()),
      xu(2, node_gd.u_h()[c_node_no].value()),          // displacements
      x_1u(2, seg_gd.u_h()[a_node_no].value()),
      x_2u(2, seg_gd.u_h()[b_node_no].value()),
      x_coord = x+xu,                                   // lagrangian(material) coordinates
      x_1_coord = x_1+x_1u,
      x_2_coord = x_2+x_2u,
      v_1 = x_coord - x_1_coord,
      v_2 = x_2_coord - x_1_coord;
   	alpha = (v_1 * v_2) / pow(norm(v_2, 2), 2);       // scalar-inner product
   if(alpha > 0.0 && alpha < 1.0) {                     // if projection onto the segment
   	C0 x_bar = (1.0-alpha)*x_1_coord + alpha*x_2_coord;
      C0 n(2,(double*)0);                                    // unit normal
      n[0] = -v_2[1]/norm(v_2,2); n[1] = v_2[0]/norm(v_2,2);
      double gap = (double)((x_coord-x_bar)*n); // projection to unit normal; inner product
   	if(gap <= tolerance_gap) {     // if within snap range(0 < gap < tolerance_gap) or penetration occured(gap < 0)
      	ret_coord = (x_bar+x_coord)/2.0;               // return contact node coordinates
      	return TRUE;
      }
   }
   return FALSE;
}
//-end of geometrical utilities-------------------------------------------------

//------------------------------------------------------------------------------
//
//                   Contact Searching Algorithm
//
//------------------------------------------------------------------------------
int setup_order_table(int *order, int size, Global_Discretization_Gamma_h_i &a_gd,
	Global_Discretization_Gamma_h_i &b_gd, double tolerance_gap) {
   Gamma_h_i &a = a_gd.gamma_h_i(), &b = b_gd.gamma_h_i();
   int a_total_node_no = a.total_element_no()+1,
       b_total_node_no = b.total_element_no()+1;
   int **order_table = new int*[2];
   order_table[0] = order;
   order_table[1] = order_table[0]+size;
   // order table example
   //  row# \ column# |    1           2          3          4          5        6       7          8       ...
   //  ----------------------------------------------------------------------------------------------------------
   //  0 (a-surface)  | a_node_#  |   -1     | a_node_# | a_node_# |    -1    | -1  |   -1     | a_node_# | ...
   //  ----------------------------------------------------------------------------------------------------------
   //  1 (b-surface)  |   -1      | b_node_# |    -1    | b_node_# | b_node_# | -1  | b_node_# |    -1    | ...
   //  ----------------------------------------------------------------------------------------------------------
   // explanation:
   //  1. first row store node # of a-surface now porjected to b-surface
   //     vice versa for the second row and a/b - surfaces;
   //  2. first tow columns show a a-node projected to b-surface, followed by
   //     a b-node porjected to a-surface
   //  3. fourth column shows an example of node-on-node projection
   //  4. sixth column with both rows as -1 indicates there is a detachment
   //     segmenet of arbitary length in the middle
   // contact node searching
   int a_node_no = 0, b_node_no = 0,
       a_node_no_cache = -1, b_node_no_cache = -1,
       order_table_count = 0;
   do {
   	Contact_Condition_Code contact_condition = contact_node_searching(a_gd, a_node_no, b_gd, b_node_no, tolerance_gap);
      // detachment segment occurred
      if(!(a_node_no_cache == -1 && b_node_no_cache == -1) &&                   // not at the beginning
         (contact_condition == A_NODE_LEADING && a_node_no > (a_node_no_cache + 1)) || // w/ any node-skipping
         (contact_condition == B_NODE_LEADING && b_node_no > (b_node_no_cache + 1))
        ) {
      	order_table[0][order_table_count] = -1;
      	order_table[1][order_table_count] = -1;
         order_table_count++;
      }                   
      a_node_no_cache = a_node_no;
      b_node_no_cache = b_node_no;
      if(contact_condition == NODE_ON_NODE) {
      	order_table[0][order_table_count] = a_node_no++;
      	order_table[1][order_table_count] = b_node_no++;
      	order_table_count++;
      } else if(contact_condition == A_NODE_LEADING) {
      	order_table[0][order_table_count] = a_node_no++;
         order_table[1][order_table_count] = -1;
      	order_table_count++;
      } else if(contact_condition == B_NODE_LEADING) {
      	order_table[0][order_table_count] = -1;
         order_table[1][order_table_count] = b_node_no++;
      	order_table_count++;
      } else if(contact_condition == PROJECTION_FAIL) {
      	// do nothing
      }
   } while(a_node_no < a_total_node_no || b_node_no < b_total_node_no);
   
   delete [] order_table;
   return order_table_count;