contact_searching.cpp

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

}

void define_contact_node(Contact_Omega_h_i &interface_oh, int *order, int size,
	int order_table_count, 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_element_no = a.total_element_no(),
       b_total_element_no = b.total_element_no();
   int **order_table = new int*[2];
   order_table[0] = order;
   order_table[1] = order_table[0]+size;
   Dynamic_Array<Node> &na = interface_oh.node_array();
   na.~Dynamic_Array<Node>(); // explicit call of the destructor
   for(int i = 0; i < order_table_count; i++) {
   	double v[2];
      Node *node;
      if(order_table[0][i] == -1 && order_table[1][i] == -1) { // do nothing skip node
      } else if(order_table[0][i] > -1 && order_table[1][i] > -1) { // node-on-node;
      	int a_node_no, b_node_no;
         if(order_table[0][i] == a_total_element_no) a_node_no = a(order_table[0][i]-1)[1];
         else a_node_no = a(order_table[0][i])[0];
         if(order_table[1][i] == b_total_element_no) b_node_no = b(order_table[1][i]-1)[1];
         else b_node_no = b(order_table[1][i])[0];
      	Node &a_node = a[a_node_no],
              &b_node = b[b_node_no];
         C0 a_u(2, a_gd.u_h()[a_node_no].value()),
            b_u(2, b_gd.u_h()[b_node_no].value());
      	v[0] = ((a_node[0]+a_u[0])+(b_node[0]+b_u[0]))/2.0; v[1] = ((a_node[1]+a_u[1])+(b_node[1]+b_u[1]))/2.0;
         node = new Node(i, 2, v);
         interface_oh.node_array().add(node);
      } else if(order_table[0][i] > -1 && order_table[1][i] == -1) { // a-node projection
      	int a_node_no;
         if(order_table[0][i] == a_total_element_no) a_node_no = a(order_table[0][i]-1)[1];
         else a_node_no = a(order_table[0][i])[0];
      	Node &a_node = a[a_node_no];
         C0 ret_coord(2, (double*)0);
         double alpha;
      	node_projection(a_node, a_gd, b_gd, alpha, tolerance_gap, 0, ret_coord);
         v[0] = (double)ret_coord[0]; v[1] = (double)ret_coord[1];
         node = new Node(i, 2, v);
         interface_oh.node_array().add(node);
      } else if(order_table[0][i] == -1 && order_table[1][i] > -1) { // b-node projection
      	int b_node_no;
         if(order_table[1][i] == b_total_element_no) b_node_no = b(order_table[1][i]-1)[1];
         else b_node_no = b(order_table[1][i])[0];
      	Node &b_node = b[b_node_no];
         C0 ret_coord(2, (double*)0);
         double alpha;
      	node_projection(b_node, b_gd, a_gd, alpha, tolerance_gap, 0, ret_coord);
         v[0] = (double)ret_coord[0]; v[1] = (double)ret_coord[1];
         node = new Node(i, 2, v);
         interface_oh.node_array().add(node);
      }
   }
   delete [] order_table;
}

void define_contact_element(Contact_Omega_h_i &interface_oh, int *order, int size,
	int order_table_count, 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_element_no = a.total_element_no(),
       b_total_element_no = b.total_element_no();
   int **order_table = new int*[2];
   order_table[0] = order;
   order_table[1] = order_table[0]+size;
	int ena[2];   // element node number array
   int aea[4];   // associate element array
   double pp[4]; // projection parameters
   double nv[4]; // normal vectors
   Omega_eh *elem;
   Contact_Omega_eh *c_element;
   int skip_node_count = 0, contact_element_count = 0;
   Dynamic_Array<Omega_eh> &ea = interface_oh.omega_eh_array();
   ea.~Dynamic_Array<Omega_eh>(); // explicit call of the destructor
   for(int i = 0; i < order_table_count-1; i++) {
      ena[0] = i - skip_node_count; ena[1] = (i - skip_node_count) + 1;
      elem = new Omega_eh(contact_element_count++, 0, 0, 2, ena);
   	if(order_table[0][i] > -1 && order_table[1][i] > -1) { // node-on-node;
         aea[0] = order_table[0][i]; aea[1] = order_table[1][i];
         pp[0] = 0.0; pp[2] = 0.0;
   		if(order_table[0][i] > -1 && order_table[1][i] > -1) {             // follow by node-on-node;
      		pp[1] = 0.0; pp[3] = 0.0;
            Node &node0 = interface_oh[contact_element_count-1], &node1 = interface_oh[contact_element_count];
            double x = node1[0]-node0[0],
                   y = node1[1]-node0[1],
                   len = sqrt(pow(x,2)+pow(y,2));
      		nv[0] = -y/len; nv[1] = x/len; nv[2] = -y/len; nv[3] = x/len; // normal vector
      		c_element = new Contact_Omega_eh(aea, pp, nv, *elem);
   			interface_oh.omega_eh_array().add(c_element);
   		} else if(order_table[0][i+1] > -1 && order_table[1][i+1] == -1) { // follow by a-node projection 
      		int a_node_no;
         	if(order_table[0][i+1] == a_total_element_no) a_node_no = a(order_table[0][i+1]-1)[1];
         	else a_node_no = a(order_table[0][i+1])[0];
      		Node &a_node = a[a_node_no];
            double alpha;
      		node_projection(a_node, a_gd, b_gd, alpha, tolerance_gap, 0, C0(2,(double*)0));
   			pp[1] = 0.0; pp[3] = 1.0-alpha;
            Node &node0 = b[b(aea[1])[0]], &node1 = b[b(aea[1])[1]];
         	int node0_no = node0.node_no(),
                node1_no = node1.node_no();
            C0 node0_u(2, b_gd.u_h()[node0_no].value()),
               node1_u(2, b_gd.u_h()[node1_no].value());
            double x = (node0[0]+node0_u[0])-(node1[0]+node1_u[0]),
                   y = (node0[1]+node0_u[1])-(node1[1]+node1_u[1]),
                   len = sqrt(pow(x,2)+pow(y,2));
      		nv[0] = -y/len; nv[1] = x/len; nv[2] = -y/len; nv[3] = x/len; // normal vector
      		c_element = new Contact_Omega_eh(aea, pp, nv, *elem);
   			interface_oh.omega_eh_array().add(c_element);
      	} else if(order_table[0][i+1] == -1 && order_table[1][i+1] > -1) { // follow by b-node projection
      		int b_node_no;
         	if(order_table[1][i+1] == b_total_element_no) b_node_no = b(order_table[1][i+1]-1)[1];
         	else b_node_no = b(order_table[1][i+1])[0];
      		Node &b_node = b[b_node_no];
            double alpha;
      		node_projection(b_node, b_gd, a_gd, alpha, tolerance_gap, 0, C0(2,(double*)0));
   			pp[1] = 1.0-alpha; pp[3] = 0.0;
            Node &node0 = a[a(aea[0])[0]], &node1 = a[a(aea[0])[1]];
         	int node0_no = node0.node_no(),
                node1_no = node1.node_no();
            C0 node0_u(2, a_gd.u_h()[node0_no].value()),
               node1_u(2, a_gd.u_h()[node1_no].value());
            double x = (node1[0]+node1_u[0])-(node0[0]+node0_u[0]),
                   y = (node1[1]+node1_u[1])-(node0[1]+node0_u[1]),
                   len = sqrt(pow(x,2)+pow(y,2));
      		nv[0] = -y/len; nv[1] = x/len; nv[2] = -y/len; nv[3] = x/len; // normal vector
      		c_element = new Contact_Omega_eh(aea, pp, nv, *elem);
   			interface_oh.omega_eh_array().add(c_element);
      	}
   	} else if(order_table[0][i] > -1 && order_table[1][i] == -1) { // a-node projection
         aea[0] = order_table[0][i];
         Node &a_node = a[a(order_table[0][i])[0]];
         double alpha;
         aea[1] = node_projection(a_node, a_gd, b_gd, alpha, tolerance_gap, 0, C0(2,(double*)0));
         pp[0] = 0.0; pp[2] = alpha;
   		if(order_table[0][i] > -1 && order_table[1][i] > -1) {             // follow by node-on-node;
         	pp[1] = pp[3] = 0.0;
            Node &node0 = b[b(aea[1])[0]], &node1 = b[b(aea[1])[1]];
         	int node0_no = node0.node_no(),
                node1_no = node1.node_no();
            C0 node0_u(2, b_gd.u_h()[node0_no].value()),
               node1_u(2, b_gd.u_h()[node1_no].value());
            double x = (node0[0]+node0_u[0])-(node1[0]+node1_u[0]),
                   y = (node0[1]+node0_u[1])-(node1[1]+node1_u[1]),
                   len = sqrt(pow(x,2)+pow(y,2));
      		nv[0] = -y/len; nv[1] = x/len; nv[2] = -y/len; nv[3] = x/len; // normal vector
      		c_element = new Contact_Omega_eh(aea, pp, nv, *elem);
   			interface_oh.omega_eh_array().add(c_element);
   		} else if(order_table[0][i+1] > -1 && order_table[1][i+1] == -1) { // follow by a-node projection
      		int a_node_no;
         	if(order_table[0][i+1] == a_total_element_no) a_node_no = a(order_table[0][i+1]-1)[1];
         	else a_node_no = a(order_table[0][i+1])[0];
      		Node &a_node = a[a_node_no];
            double alpha;
      		node_projection(a_node, a_gd, b_gd, alpha, tolerance_gap, 0, C0(2,(double*)0));
         	pp[1] = 0.0; pp[3] = 1.0-alpha;
            Node &node0 = b[b(aea[1])[0]], &node1 = b[b(aea[1])[1]];
         	int node0_no = node0.node_no(),
                node1_no = node1.node_no();
            C0 node0_u(2, b_gd.u_h()[node0_no].value()),
               node1_u(2, b_gd.u_h()[node1_no].value());
            double x = (node0[0]+node0_u[0])-(node1[0]+node1_u[0]),
                   y = (node0[1]+node0_u[1])-(node1[1]+node1_u[1]),
                   len = sqrt(pow(x,2)+pow(y,2));
      		nv[0] = -y/len; nv[1] = x/len; nv[2] = -y/len; nv[3] = x/len; // normal vector
      		c_element = new Contact_Omega_eh(aea, pp, nv, *elem);
   			interface_oh.omega_eh_array().add(c_element);
      	} else if(order_table[0][i+1] == -1 && order_table[1][i+1] > -1) { // follow by b-node projection
      		int b_node_no;
         	if(order_table[1][i+1] == b_total_element_no) b_node_no = b(order_table[1][i+1]-1)[1];
         	else b_node_no = b(order_table[1][i+1])[0];
      		Node &b_node = b[b_node_no];
            double alpha;
      		node_projection(b_node, b_gd, a_gd, alpha, tolerance_gap, 0, C0(2,(double*)0));
         	pp[1] = 1.0-alpha; pp[3] = 0.0;
            Node node0 = b[b(aea[1])[0]], node1 = b[b(aea[1])[1]];
         	int node0_no = node0.node_no(),
                node1_no = node1.node_no();
            C0 node0_u(2, b_gd.u_h()[node0_no].value()),
               node1_u(2, b_gd.u_h()[node1_no].value());
            double x = (node0[0]+node0_u[0])-(node1[0]+node1_u[0]),
                   y = (node0[1]+node0_u[1])-(node1[1]+node1_u[1]),
                   len = sqrt(pow(x,2)+pow(y,2));
      		nv[0] = -y/len; nv[1] = x/len;
            node0 = a[a(aea[0])[0]]; node1 = a[a(aea[0])[1]];
         	node0_no = node0.node_no();
            node1_no = node1.node_no();
            node0_u = C0(2, a_gd.u_h()[node0_no].value());
            node1_u = C0(2, a_gd.u_h()[node1_no].value());
            x = (node1[0]+node1_u[0])-(node0[0]+node0_u[0]);
            y = (node1[1]+node1_u[1])-(node0[1]+node0_u[1]);
            len = sqrt(pow(x,2)+pow(y,2));
            nv[2] = -y/len; nv[3] = x/len;
      		c_element = new Contact_Omega_eh(aea, pp, nv, *elem);
   			interface_oh.omega_eh_array().add(c_element);
      	}
   	} else if(order_table[0][i] == -1 && order_table[1][i] > -1) { // b-node projection
         aea[1] = order_table[1][i];
         Node &b_node = b[b(order_table[1][i])[0]];
         double alpha;
         aea[0] = node_projection(b_node, b_gd, a_gd, alpha, tolerance_gap, 0, C0(2,(double*)0));
         pp[0] = alpha; pp[2] = 0.0;
   		if(order_table[0][i] > -1 && order_table[1][i] > -1) {             // follow by node-on-node;
         	pp[1] = pp[3] = 0.0;
            Node &node0 = a[a(aea[0])[0]], &node1 = a[a(aea[0])[1]];
         	int node0_no = node0.node_no(),
                node1_no = node1.node_no();
            C0 node0_u(2, a_gd.u_h()[node0_no].value()),
               node1_u(2, a_gd.u_h()[node1_no].value());
            double x = (node1[0]+node1_u[0])-(node0[0]+node0_u[0]),
                   y = (node1[1]+node1_u[1])-(node0[1]+node0_u[1]),
                   len = sqrt(pow(x,2)+pow(y,2));
      		nv[0] = -y/len; nv[1] = x/len; nv[2] = -y/len; nv[3] = x/len; // normal vector
      		c_element = new Contact_Omega_eh(aea, pp, nv, *elem);
   			interface_oh.omega_eh_array().add(c_element);
   		} else if(order_table[0][i+1] > -1 && order_table[1][i+1] == -1) { // follow by a-node projection
      		int a_node_no;
         	if(order_table[0][i+1] == a_total_element_no) a_node_no = a(order_table[0][i+1]-1)[1];
         	else a_node_no = a(order_table[0][i+1])[0];
      		Node &a_node = a[a_node_no];
            double alpha;
      		node_projection(a_node, a_gd, b_gd, alpha, tolerance_gap, 0, C0(2,(double*)0));
         	pp[1] = 0.0; pp[3] = 1.0-alpha;
            Node node0 = a[a(aea[0])[0]], node1 = a[a(aea[0])[1]];
         	int node0_no = node0.node_no(),
                node1_no = node1.node_no();
            C0 node0_u(2, a_gd.u_h()[node0_no].value()),
               node1_u(2, a_gd.u_h()[node1_no].value());
            double x = (node1[0]+node1_u[0])-(node0[0]+node0_u[0]),
                   y = (node1[1]+node1_u[1])-(node0[1]+node0_u[1]),
                   len = sqrt(pow(x,2)+pow(y,2));
      		nv[0] = -y/len; nv[1] = x/len;
            node0 = b[b(aea[1])[0]]; node1 = b[b(aea[1])[1]];
         	node0_no = node0.node_no();
            node1_no = node1.node_no();
            node0_u = C0(2, b_gd.u_h()[node0_no].value());
            node1_u = C0(2, b_gd.u_h()[node1_no].value());
            x = (node0[0]+node0_u[0])-(node1[0]+node1_u[0]);
            y = (node0[1]+node0_u[1])-(node1[1]+node1_u[1]);
            len = sqrt(pow(x,2)+pow(y,2));
            nv[2] = -y/len; nv[3] = x/len;
      		c_element = new Contact_Omega_eh(aea, pp, nv, *elem);
   			interface_oh.omega_eh_array().add(c_element);
      	} else if(order_table[0][i+1] == -1 && order_table[1][i+1] > -1) { // follow by b-node projection
      		int b_node_no;
         	if(order_table[1][i+1] == b_total_element_no) b_node_no = b(order_table[1][i+1]-1)[1];
         	else b_node_no = b(order_table[1][i+1])[0];
      		Node &b_node = b[b_node_no];
            double alpha;
      		node_projection(b_node, b_gd, a_gd, alpha, tolerance_gap, 0, C0(2,(double*)0));
         	pp[1] = 1.0-alpha; pp[3] = 0.0;
            Node &node0 = a[a(aea[0])[0]], &node1 = a[a(aea[0])[1]];
         	int node0_no = node0.node_no(),
                node1_no = node1.node_no();
            C0 node0_u(2, a_gd.u_h()[node0_no].value()),
               node1_u(2, a_gd.u_h()[node1_no].value());
            double x = (node1[0]+node1_u[0])-(node0[0]+node0_u[0]),
                   y = (node1[1]+node1_u[1])-(node0[1]+node0_u[1]),
                   len = sqrt(pow(x,2)+pow(y,2));
      		nv[0] = -y/len; nv[1] = x/len; nv[2] = -y/len; nv[3] = x/len; // normal vector
      		c_element = new Contact_Omega_eh(aea, pp, nv, *elem);
   			interface_oh.omega_eh_array().add(c_element);
      	}
		} else if(order_table[0][i] == -1 && order_table[1][i] == -1) { // middle detachment
      	skip_node_count++;
   	}
   }
   delete [] order_table;
}