matrix_representation.cpp

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

		}
	}
   return the_global_tensor;
}

void Matrix_Representation::assembly(int nodal_load_flag, int element_no1,int element_no2, int element_incr) {
   // add nodal force (natural boundary conditions) from the rhs
   Omega_h &oh = the_global_discretization.omega_h();
	int total_node_no = oh.total_node_no(),  // no distinct number
	    ndf = the_global_discretization.u_h()[0].no_of_dim();
   Dynamic_Array<Nodal_Value>& the_u_h_array = the_global_discretization.u_h().u_h_array();
   if(nodal_load_flag) {
		for(int i = 0; i < total_node_no; i++) {
			Nodal_Constraint& gh = the_global_discretization.gh_on_gamma_h().gh_array()[i];
			for(int j = 0; j < ndf; j++) {
      		int equation_no = eqn[i][j];
   			if(gh(j) == gh_on_Gamma_h::Neumann && gh[j] != 0.0)
         		the_rhs->global_tensor()[equation_no] += gh[j];
      	}
   	}
   }
   if(Assembly_Switch == REACTION) // initialized when used
   	the_global_nodal_value &= C0(total_node_no, ndf, (double*)0);
   C0 count;
   if(Assembly_Switch == NODAL_STRESS || Assembly_Switch == NODAL_STRAIN) { // initialized when used
   	int stress_no;
      if(!axisymmetric_flag) stress_no = (ndf+1)*ndf/2;
      else stress_no = 4;
   	the_global_nodal_value &= C0(total_node_no, stress_no, (double*)0);
      if(!count.rep_ptr()) count &= C0(total_node_no, (double*)0);
      else count = 0.0;
   } else if(Assembly_Switch == NODAL_FLUX) { // HEAT CONDUCTION && IRROTATIONAL FLOW PROBLEMS
   	the_global_nodal_value &= C0(total_node_no, 2, (double*)0);
      if(!count.rep_ptr()) count &= C0(total_node_no, (double*)0);
      else count = 0.0;
   } else if(Assembly_Switch == NODAL_SCALAR) { // yield ratio
   	the_global_nodal_value &= C0(total_node_no, (double*)0);
      if(!count.rep_ptr()) count &= C0(total_node_no, (double*)0);
      else count = 0.0;
   }
   // add stiffness, reaction (essential boundary conditions) and distributed load
	if(element_no2 == -1) element_no2 = element_no1; // assemble only one element
	if(element_no1 == -1) {                          // assemble all distinct elements
   	int total_element_no = oh.omega_eh_array().length();
		for(int i = 0; i < total_element_no; i++) {
			Omega_eh& elem = oh(i);
         int etn = elem.element_type_no();
         Element_Formulation *element_type = Element_Formulation::type_list;
         if(etn > 0) for(int j = 0; j < etn; j++) element_type = element_type->next;
         Element_Formulation ef = element_type->create(i, the_global_discretization); // form element
         if(Assembly_Switch == ALL || Assembly_Switch == LHS) {
				// a. add element stiffness matrix
				Element_Tensor element_lhs_matrix(i, ef.lhs(), *this);
				(*the_lhs) += element_lhs_matrix;
         }
         if(Assembly_Switch == ALL || Assembly_Switch == RHS) {
				// b. add reaction and distributed load
				Element_Tensor element_rhs_vector(i, ef.rhs(), *this);
				rhs() += element_rhs_vector;
         }
         if(Assembly_Switch == REACTION) {
            int nen = elem.element_node_no();
            C0 element_reaction = ef.element_nodal_value();
            for(int i = 0; i < nen; i++) {
            	int node_no = elem[i],
                   tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
               if(tie_node_no != -1) node_no = tie_node_no;
            	for(int j = 0; j < ndf; j++)
						the_global_nodal_value[oh.node_order(node_no)][j] += element_reaction[i*ndf+j];
            }
         } else if(Assembly_Switch == NODAL_STRESS || Assembly_Switch == NODAL_STRAIN) {
            int nen = elem.element_node_no();
            C0 element_stress = ef.element_nodal_value();
   			int stress_no;
      		if(!axisymmetric_flag) stress_no = (ndf+1)*ndf/2;
      		else stress_no = 4;
            for(int i = 0; i < nen; i++) {
            	int node_no = elem[i],
                   tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
               if(tie_node_no != -1) node_no = tie_node_no;
            	for(int j = 0; j < stress_no; j++)
						the_global_nodal_value[oh.node_order(node_no)][j] += element_stress[i*stress_no+j];
               count[oh.node_order(node_no)] += 1;
            }
         } else if (Assembly_Switch == NODAL_FLUX) {
            int nen = elem.element_node_no();
            C0 element_velocity = ef.element_nodal_value();
            for(int i = 0; i < nen; i++) {
            	int node_no = elem[i],
                   tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
               if(tie_node_no != -1) node_no = tie_node_no;
            	for(int j = 0; j < 2; j++)
						the_global_nodal_value[oh.node_order(node_no)][j] += element_velocity[i*2+j];
               count[oh.node_order(node_no)] += 1;
            }
         } else if (Assembly_Switch == NODAL_SCALAR) {
            int nen = elem.element_node_no();
            C0 element_scalar = ef.element_nodal_value();
            for(int i = 0; i < nen; i++) {
            	int node_no = elem[i],
                   tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
               if(tie_node_no != -1) node_no = tie_node_no;
            	the_global_nodal_value[oh.node_order(node_no)] += element_scalar[i];
               count[oh.node_order(node_no)] += 1;
            }
         }
		}
	} else { // assemble specified range of elements
   	for(int i = element_no1; i <= element_no2; i+=element_incr) {
   		Omega_eh& elem = oh(i);
      	int etn = elem.element_type_no();
      	Element_Formulation *element_type = Element_Formulation::type_list;
      	if(etn > 0) for(int j = 0; j < etn; j++) element_type = element_type->next;
      	Element_Formulation ef = element_type->create(i, the_global_discretization); // form element
         if(Assembly_Switch == ALL || Assembly_Switch == LHS) {
				// a. add element stiffness matrix
				Element_Tensor element_lhs_matrix(i, ef.lhs(), *this);
				(*the_lhs) += element_lhs_matrix;
         }
         if(Assembly_Switch == ALL || Assembly_Switch == RHS) {
				// b. add reaction and distributed load
				Element_Tensor element_rhs_vector(i, ef.rhs(), *this);
				rhs() += element_rhs_vector;
         }
         if(Assembly_Switch == REACTION) {
            int nen = elem.element_node_no();
            C0 element_reaction = ef.element_nodal_value();
            for(int i = 0; i < nen; i++) {
            	int node_no = elem[i],
                   tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
               if(tie_node_no != -1) node_no = tie_node_no;
            	for(int j = 0; j < ndf; j++)
						the_global_nodal_value[oh.node_order(node_no)][j] += element_reaction[i*ndf+j];
            }
         } else if(Assembly_Switch == NODAL_STRESS || Assembly_Switch == NODAL_STRAIN) {
            int nen = elem.element_node_no();
            C0 element_stress = ef.element_nodal_value();
   			int stress_no;
      		if(!axisymmetric_flag) stress_no = (ndf+1)*ndf/2;
      		else stress_no = 4;
            for(int i = 0; i < nen; i++) {
            	int node_no = elem[i],
                   tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
               if(tie_node_no != -1) node_no = tie_node_no;
            	for(int j = 0; j < stress_no; j++)
						the_global_nodal_value[oh.node_order(node_no)][j] += element_stress[i*stress_no+j];
               count[oh.node_order(node_no)] += 1;
            }
         } else if (Assembly_Switch == NODAL_FLUX) {
            int nen = elem.element_node_no();
            C0 element_velocity = ef.element_nodal_value();
            for(int i = 0; i < nen; i++) {
            	int node_no = elem[i],
                   tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
               if(tie_node_no != -1) node_no = tie_node_no;
            	for(int j = 0; j < 2; j++)
						the_global_nodal_value[oh.node_order(node_no)][j] += element_velocity[i*2+j];
               count[oh.node_order(node_no)] += 1;
            }
         } else if (Assembly_Switch == NODAL_SCALAR) {
            int nen = elem.element_node_no();
            C0 element_scalar = ef.element_nodal_value();
            for(int i = 0; i < nen; i++) {
            	int node_no = elem[i],
                   tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
               if(tie_node_no != -1) node_no = tie_node_no;
            	the_global_nodal_value[oh.node_order(node_no)] += element_scalar[i];
               count[oh.node_order(node_no)] += 1;
            }
         }
	 	}
   }
   if(Assembly_Switch == NODAL_STRESS || Assembly_Switch == NODAL_STRAIN) {
   	for(int i = 0; i < total_node_no; i++) {
      	int node_no = oh.node_array()[i].node_no(),
             tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
         if(tie_node_no != -1) node_no = tie_node_no;
      	if(fabs((double)count[i]) >= 1.e-15) {
   			int stress_no;
      		if(!axisymmetric_flag) stress_no = (ndf+1)*ndf/2;
      		else stress_no = 4;
         	for(int j = 0; j < stress_no; j++)
            	the_global_nodal_value[oh.node_order(node_no)][j] /= (double)count[oh.node_order(node_no)];
         }
      }
   } else if (Assembly_Switch == NODAL_FLUX) {
   	for(int i = 0; i < total_node_no; i++) {
      	int node_no = oh.node_array()[i].node_no(),
             tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
         if(tie_node_no != -1) node_no = tie_node_no;
      	if(fabs((double)count[i]) >= 1.e-15) {
         	for(int j = 0; j < 2; j++)
            	the_global_nodal_value[oh.node_order(node_no)][j] /= (double)count[oh.node_order(node_no)];
         }
      }      
   } else if (Assembly_Switch == NODAL_SCALAR) {
   	for(int i = 0; i < total_node_no; i++) {
      	int node_no = oh.node_array()[i].node_no(),
             tie_node_no = the_u_h_array[oh.node_order(node_no)].tie_node_no();
         if(tie_node_no != -1) node_no = tie_node_no;
      	if(fabs((double)count[i]) >= 1.e-15) {
         	the_global_nodal_value[oh.node_order(node_no)] /= (double)count[oh.node_order(node_no)];
         }
      }
   }
}