block.cpp

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

#include "vs_h.h"
#include "dynamic_array.h"
#include "omega_h.h"
#include "block.h"

void block(Omega_h* oh, int row_node_no, int column_node_no, int
	control_node_no, int* control_node_flag, double* coord,
   int first_node_no, int first_element_no, int column_node_skip, int column_element_skip,
   int element_type_id, int material_type_id) {
   double v[2];
   C0 x(control_node_no, 2, coord), N(control_node_no, (double*)0);
	for(int i = 0; i < row_node_no; i++)
   for(int j = 0; j < column_node_no; j++) {
   	int nn = first_node_no + i*(column_node_no+column_node_skip) + j;
      v[0] = 2.0*(double)j/(double)(column_node_no-1)-1.0; // natural coordinates
      v[1] = 2.0*(double)i/(double)(row_node_no-1)-1.0;
      shape(v, control_node_no, control_node_flag, N);
   	C0 X = N*x;                                           // physical Coordinates
      v[0] = X[0]; v[1] = X[1];
      Node* node = new Node(nn, 2, v);
      oh->set(node);
   }
   switch(EP::ep) {
   	case EP::QUADRILATERALS_4_NODES:
      	q4(oh->omega_eh_array(), row_node_no, column_node_no, first_node_no, first_element_no,
         	column_node_skip, column_element_skip, element_type_id, material_type_id);
      	break;
      case EP::LAGRANGIAN_9_NODES:
      	l9(oh->omega_eh_array(), row_node_no, column_node_no, first_node_no, first_element_no,
         	column_node_skip, column_element_skip, element_type_id, material_type_id);
      	break;
      case EP::SLASH_TRIANGLES:
      	s_tri(oh->omega_eh_array(), row_node_no, column_node_no, first_node_no, first_element_no,
         	column_node_skip, column_element_skip, element_type_id, material_type_id);
      	break;
      case EP::BACKSLASH_TRIANGLES:
      	b_tri(oh->omega_eh_array(), row_node_no, column_node_no, first_node_no, first_element_no,
         	column_node_skip, column_element_skip, element_type_id, material_type_id);
      	break;
      case EP::SB_ALTERNATE_TRIANGLES:
      	sb_tri(oh->omega_eh_array(), row_node_no, column_node_no, first_node_no, first_element_no,
         	column_node_skip, column_element_skip, element_type_id, material_type_id);
      	break;
      case EP::BS_ALTERNATE_TRIANGLES:
      	bs_tri(oh->omega_eh_array(), row_node_no, column_node_no, first_node_no, first_element_no,
         	column_node_skip, column_element_skip, element_type_id, material_type_id);
      	break;
   }
}

void shape(double* v, int control_node_no, int* control_node_flag, C0& N) {
	C0 Zai(v[0]), Eta(v[1]);
   // shape function of 4 to 9 nodes quadrilateral
   // initial four corner nodes
   N[0] = (1.-Zai)*(1.-Eta)/4.; N[1] = (1.+Zai)*(1.-Eta)/4.;
   N[2] = (1.+Zai)*(1.+Eta)/4.; N[3] = (1.-Zai)*(1.+Eta)/4.;
   if(control_node_no == 9 && control_node_flag[8]) {
      // add center node
      N[8] = (1.-Zai.pow(2))*(1.-Eta.pow(2));
      // modification to four corner nodes due to the presence of the center node
      N[0] -= N[8]/4.; N[1] -= N[8]/4.; N[2] -= N[8]/4.; N[3] -= N[8]/4.;
   }
   if(control_node_no > 4) {
      if(control_node_no >= 5) {
      	if(control_node_flag[4]) {
         	// add edge node
            N[4] = (1.-Zai.pow(2))*(1.-Eta)/2.;
            // modification of the edge node due to the presence of center node
            if(control_node_no == 9) N[4] -= N[8]/2.;
            // modification of corner nodes due to the presence of the edge node
            N[0] -= N[4] / 2.; N[1] -= N[4]/2.;
         } else N[4] = 0.0;
      }
      if(control_node_no >= 6) {
         if(control_node_flag[5]) {
            // add edge node
            N[5] = (1.-Eta.pow(2))*(1.+Zai)/2.;
            // modification of the edge node due to the presence of center node
            if(control_node_no == 9) N[5] -= N[8]/2.;
            // modification of corner nodes due to the presence of the edge node
            N[1] -= N[5] / 2.; N[2] -= N[5]/2.;
         } else N[5] = 0.0;
      }
      if(control_node_no >= 7) {
         if(control_node_flag[6]) {
            // add edge node
            N[6] = (1.-Zai.pow(2))*(1.+Eta)/2.;
            // modification of the edge node due to the presence of center node
            if(control_node_no == 9) N[6] -= N[8]/2.;
            // modification of corner nodes due to the presence of the edge node
            N[2] -= N[6] / 2.; N[3] -= N[6]/2.;
         } else N[6] = 0.0;
      }
      if(control_node_no >= 8) {
         if(control_node_flag[7]) {
            // add edge node
            N[7] = (1.-Eta.pow(2))*(1.-Zai)/2.;
            // modification of the edge node due to the presence of center node
            if(control_node_no == 9) N[7] -= N[8]/2.;
            // modification of corner nodes due to the presence of the edge node
            N[3] -= N[7] / 2.; N[0] -= N[7]/2.;
         } else N[7] = 0.0;
      }
   }

}

void q4(Dynamic_Array<Omega_eh>& element_array, int row_node_no, int column_node_no,
	int first_node_no, int first_element_no, int column_node_skip, int column_element_skip,
   int element_type_id, int material_type_id) {
   int row_element_no = row_node_no -1,
       column_element_no = column_node_no -1;
   for(int i = 0; i < row_element_no; i++)
   for(int j = 0; j < column_element_no; j++) {
       int nn = first_node_no + i*(column_node_no+column_node_skip) + j; // node number at lower left corner
       int ena[4]; ena[0] = nn; ena[1] = ena[0]+1; ena[3] = nn + column_node_no+column_node_skip; ena[2] = ena[3]+1;
       int en = first_element_no + i*(column_element_no+column_element_skip)+j;
       Omega_eh* elem = new Omega_eh(en, element_type_id, material_type_id, 4, ena);
       element_array.add(elem);
   }
}

void l9(Dynamic_Array<Omega_eh>& element_array, int row_node_no, int column_node_no,
	int first_node_no, int first_element_no, int column_node_skip, int column_element_skip,
   int element_type_id, int material_type_id) {
   // define elements
   int row_element_no = (row_node_no-1)/2,
       column_element_no = (column_node_no-1)/2;
   for(int i = 0; i < row_element_no; i++)
   for(int j = 0; j < column_element_no; j++) {
       int nn = first_node_no + 2*i*(column_node_no+column_node_skip) + 2*j; // node number at lower left corner
       int ena[9];
       ena[0] = nn; ena[1] = ena[0]+2; ena[3] = nn + 2*(column_node_no+column_node_skip); ena[2] = ena[3]+2;
       ena[4] = ena[0] + 1; ena[5] = ena[1]+column_node_no+column_node_skip; ena[6] = ena[3]+1;
       ena[7] = ena[0] + column_node_no+column_node_skip; ena[8] = ena[7]+1;
       int en = first_element_no + i*(column_element_no+column_element_skip)+j;
       Omega_eh* elem = new Omega_eh(en, element_type_id, material_type_id, 9, ena);
       element_array.add(elem);
   }
}

void s_tri(Dynamic_Array<Omega_eh>& element_array, int row_node_no, int column_node_no,
	int first_node_no, int first_element_no, int column_node_skip, int column_element_skip,
   int element_type_id, int material_type_id) {
   int row_segment_no = row_node_no -1,
       column_segment_no = column_node_no -1;
   for(int i = 0; i < row_segment_no; i++)
   for(int j = 0; j < column_segment_no; j++) {
       int nn = first_node_no + i*(column_node_no+column_node_skip) + j; // node number at lower left corner
       int ena[3]; ena[0] = nn; ena[1] = nn + column_node_no+column_node_skip+1; ena[2] = ena[1]-1;
       int en = first_element_no + i*(column_segment_no*2+column_element_skip)+j*2;
       Omega_eh* elem = new Omega_eh(en, element_type_id, material_type_id, 3, ena);
       element_array.add(elem);
       ena[0] = nn; ena[1] = ena[0]+1; ena[2] = nn + column_node_no+column_node_skip+1;
       elem = new Omega_eh(en+1, element_type_id, material_type_id, 3, ena);
       element_array.add(elem);
   }
}

void b_tri(Dynamic_Array<Omega_eh>& element_array, int row_node_no, int column_node_no,
	int first_node_no, int first_element_no, int column_node_skip, int column_element_skip,
   int element_type_id, int material_type_id) {
   int row_segment_no = row_node_no -1,
       column_segment_no = column_node_no -1;
   for(int i = 0; i < row_segment_no; i++)
   for(int j = 0; j < column_segment_no; j++) {
       int nn = first_node_no + i*(column_node_no+column_node_skip) + j; // node number at lower left corner
       int ena[3]; ena[0] = nn; ena[1] = nn+1; ena[2] = nn + column_node_no+column_node_skip;
       int en = first_element_no + i*(column_segment_no*2+column_element_skip)+j*2;
       Omega_eh* elem = new Omega_eh(en, element_type_id, material_type_id, 3, ena);
       element_array.add(elem);
       ena[0] = nn+1; ena[1] = nn + column_node_no+column_node_skip+1; ena[2] = ena[1] -1;
       elem = new Omega_eh(en+1, element_type_id, material_type_id, 3, ena);
       element_array.add(elem);
   }
}

void sb_tri(Dynamic_Array<Omega_eh>& element_array, int row_node_no, int column_node_no,
	int first_node_no, int first_element_no, int column_node_skip, int column_element_skip,
   int element_type_id, int material_type_id) {
   int row_segment_no = row_node_no -1,
       column_segment_no = column_node_no -1;
   double iptr;
   for(int i = 0; i < row_segment_no; i++)
   for(int j = 0; j < column_segment_no; j++) {
   	if( fabs( modf(((double)j/2.0), &iptr) ) < 1.e-12) { // even number slashed triangles
      	int nn = first_node_no + i*(column_node_no+column_node_skip) + j; // node number at lower left corner
       	int ena[3]; ena[0] = nn; ena[1] = nn + column_node_no+column_node_skip+1; ena[2] = ena[1]-1;
       	int en = first_element_no + i*(column_segment_no*2+column_element_skip)+j*2;
       	Omega_eh* elem = new Omega_eh(en, element_type_id, material_type_id, 3, ena);
       	element_array.add(elem);
       	ena[0] = nn; ena[1] = ena[0]+1; ena[2] = nn + column_node_no+column_node_skip+1;
       	elem = new Omega_eh(en+1, element_type_id, material_type_id, 3, ena);
       	element_array.add(elem);
      } else { // odd number back-slashed triangles
       	int nn = first_node_no + i*(column_node_no+column_node_skip) + j; // node number at lower left corner
       	int ena[3]; ena[0] = nn; ena[1] = nn+1; ena[2] = nn + column_node_no+column_node_skip;
       	int en = first_element_no + i*(column_segment_no*2+column_element_skip)+j*2;
       	Omega_eh* elem = new Omega_eh(en, element_type_id, material_type_id, 3, ena);
       	element_array.add(elem);
       	ena[0] = nn+1; ena[1] = nn + column_node_no+column_node_skip+1; ena[2] = ena[1] -1;
       	elem = new Omega_eh(en+1, element_type_id, material_type_id, 3, ena);
       	element_array.add(elem);
      }
   }
}

void bs_tri(Dynamic_Array<Omega_eh>& element_array, int row_node_no, int column_node_no,
	int first_node_no, int first_element_no, int column_node_skip, int column_element_skip,
   int element_type_id, int material_type_id) {
   int row_segment_no = row_node_no -1,
       column_segment_no = column_node_no -1;
   double iptr;
   for(int i = 0; i < row_segment_no; i++)
   for(int j = 0; j < column_segment_no; j++) {
   	if( fabs( modf(((double)j/2.0), &iptr) ) < 1.e-12) { // even number back-slashed triangles
       	int nn = first_node_no + i*(column_node_no+column_node_skip) + j; // node number at lower left corner
       	int ena[3]; ena[0] = nn; ena[1] = nn+1; ena[2] = nn + column_node_no+column_node_skip;
       	int en = first_element_no + i*(column_segment_no*2+column_element_skip)+j*2;
       	Omega_eh* elem = new Omega_eh(en, element_type_id, material_type_id, 3, ena);
       	element_array.add(elem);
       	ena[0] = nn+1; ena[1] = nn + column_node_no+column_node_skip+1; ena[2] = ena[1] -1;
       	elem = new Omega_eh(en+1, element_type_id, material_type_id, 3, ena);
       	element_array.add(elem);
      } else { // odd number slashed triangles
      	int nn = first_node_no + i*(column_node_no+column_node_skip) + j; // node number at lower left corner
       	int ena[3]; ena[0] = nn; ena[1] = nn + column_node_no+column_node_skip+1; ena[2] = ena[1]-1;
       	int en = first_element_no + i*(column_segment_no*2+column_element_skip)+j*2;
       	Omega_eh* elem = new Omega_eh(en, element_type_id, material_type_id, 3, ena);
       	element_array.add(elem);
       	ena[0] = nn; ena[1] = ena[0]+1; ena[2] = nn + column_node_no+column_node_skip+1;
       	elem = new Omega_eh(en+1, element_type_id, material_type_id, 3, ena);
       	element_array.add(elem);
      }
   }
}