omega_h_n.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 "u_h.h"
#include "element_formulation.h"
#include "matrix_representation.h"
#include "omega_h_n.h"

Global_Discretization_Couple& Global_Discretization_Couple::operator=(const Global_Discretization_Couple& a) {
	if(this == &a) return *this;
   omega_h() = a.principal_discretization().omega_h();                                             // base class members
   u_h() = a.principal_discretization().u_h();
   gh_on_gamma_h() = a.principal_discretization().gh_on_gamma_h();
   the_subordinate_discretization = a.the_subordinate_discretization; // derived class member
   type_id = a.type_id;
   return *this;
}

static Global_Discretization default_gd(Global_Discretization());
Element_Formulation_Couple::Element_Formulation_Couple(Element_Type_Register a)
   	: Element_Formulation(a), the_subordinate_discretization(default_gd) {}
// for diagonal matrix element formulation
Element_Formulation_Couple::Element_Formulation_Couple(int en, Global_Discretization& gd) :
	Element_Formulation(en, gd), the_subordinate_discretization(gd) { }

// for off-diagonal matrix element formulation
Element_Formulation_Couple::Element_Formulation_Couple(int en,
	Global_Discretization_Couple& gdc)
   : Element_Formulation(en, gdc.principal_discretization()),
   s_element_no(en), the_subordinate_discretization(gdc.subordinate_discretization()) {
   Omega_eh& s_element = the_subordinate_discretization.omega_h()(en);
   s_nen = s_element.element_node_no(),
   s_nsd = the_subordinate_discretization.omega_h()[ s_element[0] ].no_of_dim();
   s_ndf = the_subordinate_discretization.u_h()[s_element[0]].no_of_dim();
   s_material_type_no = s_element.material_type_no();
   s_xl &= the_subordinate_discretization.element_coordinate(en);
   s_gl &= the_subordinate_discretization.element_fixed_variable(en);
   s_ul &= the_subordinate_discretization.element_free_variable(en);
}

// protected virtual member function
C0& Element_Formulation_Couple::__lhs() {
	the_lhs &= stiff;
   return the_lhs;
}

//==============================================================================
// RHS definition for diagonal and off-diagonal
//
// A. diagonal-block
// call through ef->rhs() in Element_Formulation class =>
// and virtual mechanism of __rhs() to forwarding to __rhs() in the
// Element_Formulation_Couple
C0& Element_Formulation_Couple::__reaction() {
	the_reaction &= -stiff * gl;
   return the_reaction;
}

C0& Element_Formulation_Couple::__rhs() {
	the_rhs &= __reaction();
	if(force.rep_ptr()) the_rhs += force;
   return the_rhs;
}

// B. off-diagonal
// call through explicit casting ((Element_Formulation_Couple*)ef)->rhs()
// to intercept rhs() through Element_Formulation class with a different
// rhs implementation
C0& Element_Formulation_Couple::rhs() {
	the_rhs &= -stiff * s_gl;
	if(force.rep_ptr()) the_rhs += force;
   return the_rhs;
}
//==============================================================================

C0& Element_Formulation_Couple::__s_reaction() {
	the_s_reaction &= -(~stiff) * gl;
   return the_s_reaction;
}

C0& Element_Formulation_Couple::__s_rhs() {
	the_s_rhs &= __s_reaction();
	if(s_force.rep_ptr()) the_s_rhs += s_force;
   return the_s_rhs;
}

void Matrix_Representation_Couple::__initialization(char *s, Global_Tensor *principal_rhs, Global_Tensor *subordinate_rhs) {
	// setup subordinate dof variable table
	int total_node_no = the_subordinate_discretization.omega_h().total_node_no(), // no distinct number
	    ndf = the_subordinate_discretization.u_h()[0].no_of_dim();
	var = new int*[total_node_no];
	var[0] = new int[total_node_no*ndf];
	for(int i = 1; i < total_node_no; i++) var[i] = var[i-1] + ndf;
	int count = 0;
	for(i = 0; i < total_node_no; i++)
	for(int j = 0; j < ndf; j++) {
		Nodal_Constraint& gh = the_subordinate_discretization.gh_on_gamma_h()[i];
		if(gh(j) == gh_on_Gamma_h::Dirichlet) var[i][j] = Matrix_Representation::Fixed_Variable; // fixed degree of freedom
		else if(gh.tie_node_no() != -1) var[i][j] = Matrix_Representation::Tied_Node;            // tie node
		else var[i][j] = count++;
	}
   total_var_no = count;
   rhs_owner = FALSE; // initial state
   if(!s) { // default use non-sparse matrix for rapid iteration of a small size problem
      the_lhs = new Coupled_Global_Tensor(C0(total_eqn_no, total_var_no, (double*)0), *this);
      if(!principal_rhs) {
      	rhs_owner = TRUE;
      	the_rhs = new Global_Tensor(C0(total_eqn_no, (double*)0), *this);
      } else {
      	the_rhs = principal_rhs;
      }
   } else { the_lhs = 0; the_rhs = 0; } // the initial state
   the_subordinate_rhs = subordinate_rhs; // always referenced to an existing rhs
}

Matrix_Representation_Couple::Matrix_Representation_Couple(
	Global_Discretization_Couple& gdc, char *s,
   Global_Tensor* principal_rhs, Global_Tensor* subordinate_rhs,
   Matrix_Representation* subordinate_representation) :
   Matrix_Representation(gdc.principal_discretization(), "NONE_LHS_RHS"), // pass an invalid string to avoid forming lhs & rhs
   the_subordinate_discretization(gdc.subordinate_discretization()),
   gdc_type_id((Global_Discretization_Couple*)(gdc.type())),
   the_subordinate_matrix_representation(subordinate_representation) {
   __initialization(s, principal_rhs, subordinate_rhs);
}

Matrix_Representation_Couple::~Matrix_Representation_Couple() {
	if(var[0] && var) {
   	delete [] var[0];
      var[0] = 0;
   }
   if(var) {
   	delete [] var;
      var = 0;
   }
   if(the_lhs) {
   	delete the_lhs;
      the_lhs = 0;
   }
   if(the_rhs && rhs_owner) {
   	delete the_rhs;
      the_rhs = 0;
      rhs_owner = FALSE;
   }
}

Matrix_Representation_Couple::Matrix_Representation_Couple(const Matrix_Representation_Couple& a)
	: the_subordinate_discretization(a.the_subordinate_discretization), total_var_no(a.total_var_no) {
   the_global_discretization = a.the_global_discretization;
   total_eqn_no = a.total_eqn_no;
  	int total_node_no = the_global_discretization.omega_h().total_node_no(), // no distinct number
	    ndf = the_global_discretization.u_h()[0].no_of_dim();
	eqn = new int*[total_node_no];
	eqn[0] = a.eqn[0];
	for(int i = 1; i < total_node_no; i++) eqn[i] = eqn[i-1] + ndf;
  	int total_s_node_no = the_subordinate_discretization.omega_h().total_node_no(), // no distinct number
	    s_ndf = the_subordinate_discretization.u_h()[0].no_of_dim();
	var = new int*[total_s_node_no];
	var[0] = a.var[0];
	for(i = 1; i < total_s_node_no; i++) var[i] = var[i-1] + s_ndf;
   *(the_lhs) = *(a.the_lhs);
   *(the_rhs) = *(a.the_rhs);
   *(the_subordinate_rhs) = *(a.the_subordinate_rhs);
}

Matrix_Representation_Couple& Matrix_Representation_Couple::operator=(const Matrix_Representation_Couple& a) {
	if(this == &a) return *this;
	global_discretization() = a.global_discretization();
   the_subordinate_discretization = a.the_subordinate_discretization;
   total_eqn_no = a.total_eqn_no;
   int total_node_no = a.global_discretization().omega_h().total_node_no(), // no distinct number
	    ndf = a.global_discretization().u_h()[0].no_of_dim();
   if(eqn[0] && eqn) { delete [] eqn[0]; eqn[0] = 0; }
   if(eqn) { delete [] eqn; eqn = 0; }
   eqn = new int*[total_node_no];
   eqn[0] = a.eqn[0];
   for(int i = 1; i < total_node_no; i++) eqn[i] = eqn[i-1] + ndf;
   total_var_no = a.total_var_no;
   int total_s_node_no = a.the_subordinate_discretization.omega_h().total_node_no(), // no distinct number
	    s_ndf = a.the_subordinate_discretization.u_h()[0].no_of_dim();
   if(var[0] && var) { delete [] var[0]; var[0] = 0; }
   if(var) { delete [] var; var = 0; }
   var = new int*[total_s_node_no];
   var[0] = a.var[0];
   for(i = 1; i < total_s_node_no; i++) var[i] = var[i-1] + s_ndf;
   *(the_lhs) = *(a.the_lhs);
   *(the_rhs) = *(a.the_rhs);
   *(the_subordinate_rhs) = *(a.the_subordinate_rhs);
   return *this;
}
				
void Matrix_Representation_Couple::assembly(int element_no1,int element_no2, int element_incr, int nodal_load_flag) {
	// add nodal force (natural boundary conditions) from the principal rhs
	int total_node_no = (Matrix_Representation::global_discretization()).omega_h().total_node_no(), // no distinct number
	    ndf = (Matrix_Representation::global_discretization()).u_h()[0].no_of_dim();
   if(nodal_load_flag) {
		for(int i = 0; i < total_node_no; i++) {
			Nodal_Constraint& gh = (Matrix_Representation::global_discretization()).gh_on_gamma_h()[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];
      	}
   	}
   }
   // 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 = principal_discretization().omega_h().omega_eh_array().length();
		for(int i = 0; i < total_element_no; i++) {
			Omega_eh& elem = principal_discretization().omega_h()(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;
         ef.rep() = ((Element_Formulation_Couple*)element_type)->make(i,        // form element
         	Global_Discretization_Couple(principal_discretization(),
                                         the_subordinate_discretization,
                                         gdc_type_id)
            );
         if(Assembly_Switch == LHS || Assembly_Switch == LHS_AND_RHS || Assembly_Switch == ALL) {
				// a. add element stiffness matrix
				Element_Tensor element_lhs_matrix(i, ((Element_Formulation_Couple*)(ef.rep()))->lhs(), *this);
				(*the_lhs) += element_lhs_matrix;
         }
         if(Assembly_Switch == RHS || Assembly_Switch == LHS_AND_RHS || Assembly_Switch == ALL) {
				// b. add reaction and distributed load
				Element_Tensor element_rhs_vector(i, ((Element_Formulation_Couple*)(ef.rep()))->rhs(), *this);
				(*the_rhs) += element_rhs_vector;
         }
         if(Assembly_Switch == SUBORDINATE_RHS || Assembly_Switch == ALL) {
				// c. add reaction and distributed load
            Matrix_Representation *smr;
            if(!subordinate_matrix_representation()) smr = this;
            else smr = subordinate_matrix_representation();
				Element_Tensor element_subordinate_rhs_vector(i, ((Element_Formulation_Couple*)(ef.rep()))->subordinate_rhs(),
            										*smr);
            (*the_subordinate_rhs) += element_subordinate_rhs_vector;
         }
		}
	} else { // assemble specified range of elements
   	for(int i = element_no1; i <= element_no2; i+=element_incr) {
			Omega_eh& elem = principal_discretization().omega_h()(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;
         ef.rep() = ((Element_Formulation_Couple*)element_type)->make(i,        // form element
         	Global_Discretization_Couple(principal_discretization(),
                                         the_subordinate_discretization,
                                         gdc_type_id)
            );
         if(Assembly_Switch == LHS || Assembly_Switch == LHS_AND_RHS || Assembly_Switch == ALL) {
				// a. add element stiffness matrix
				Element_Tensor element_lhs_matrix(i, ((Element_Formulation_Couple*)(ef.rep()))->lhs(), *this);
				(*the_lhs) += element_lhs_matrix;
         }
         if(Assembly_Switch == RHS || Assembly_Switch == LHS_AND_RHS || Assembly_Switch == ALL) {
				// b. add reaction and distributed load
				Element_Tensor element_rhs_vector(i, ((Element_Formulation_Couple*)(ef.rep()))->rhs(), *this);
				(*the_rhs) += element_rhs_vector;
         }
         if(Assembly_Switch == SUBORDINATE_RHS || Assembly_Switch == ALL) {
				// c. add reaction and distributed load
				Element_Tensor element_subordinate_rhs_vector(i, ((Element_Formulation_Couple*)(ef.rep()))->subordinate_rhs(), *this);
				(*the_subordinate_rhs) += element_subordinate_rhs_vector;
         }
	 	}
   }
}

C0& Coupled_Global_Tensor::operator+=(Element_Tensor& element_tensor) {
   int en = element_tensor.element_no();
   Omega_eh elem = the_matrix_representation_couple.global_discretization().omega_h()(en);
	int nen = elem.element_node_no(),
       first_node = elem[0],
       ndf = the_matrix_representation_couple.global_discretization().u_h()[first_node].no_of_dim();
   if(element_tensor.element_tensor().type() == vector_type) {
   	for(int i = 0; i < nen; i++) {
         int node_no = elem[i];
         for(int j = 0; j < ndf; j++) {
            int eq_no = the_matrix_representation_couple.eqn[node_no][j];
            if(eq_no == Matrix_Representation::Tied_Node) {
               int tie_node_no = the_matrix_representation_couple.global_discretization().u_h().u_h_array()[node_no].tie_node_no();
               eq_no = the_matrix_representation_couple.eqn[tie_node_no][j];
            }
            if(eq_no >=0 )
            global_tensor()[eq_no] += element_tensor.element_tensor()[i*ndf+j];
         }
   	}
   } else if(element_tensor.element_tensor().type() == matrix_type) {
      Omega_eh s_elem = ((Matrix_Representation_Couple&)(element_tensor.matrix_representation())).
      						subordinate_discretization().omega_h()(en);
      int s_nen = s_elem.element_node_no(),
          s_first_node = s_elem[0],
          s_ndf = the_matrix_representation_couple.the_subordinate_discretization.u_h()[s_first_node].no_of_dim();
   	for(int i = 0; i < nen; i++)
      for(int j = 0; j < s_nen; j++) {
         int node_no1 = elem[i], node_no2 = s_elem[j];
         for(int k = 0; k < ndf; k++)
         for(int l = 0; l < s_ndf; l++) {
            int eqn1 = the_matrix_representation_couple.eqn[node_no1][k],
                eqn2 = the_matrix_representation_couple.var[node_no2][l];
            // converts tied node
            if(eqn1 == Matrix_Representation::Tied_Node) {
               int tie_node_no = the_matrix_representation_couple.global_discretization().u_h().u_h_array()[node_no1].tie_node_no();
               eqn1 = the_matrix_representation_couple.eqn[tie_node_no][k];
            }
            if(eqn2 == Matrix_Representation::Tied_Node) {
            	int tie_node_no = the_matrix_representation_couple.the_subordinate_discretization.u_h().u_h_array()[node_no2].tie_node_no();
               eqn2 = the_matrix_representation_couple.var[tie_node_no][l];
            }
            if(eqn1 >= 0 && eqn2 >= 0) // if none of the two dof is a fixed degree of freedom
					global_tensor()[eqn1][eqn2] += element_tensor.element_tensor()[i*ndf+k][j*s_ndf+l];
			}
		}
	}
   return global_tensor();
}