matrix_representation.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"

// private utility for constructor initialization
// form a clique of unfortunate dependency between classes U_h and Matrix_Representation
void Matrix_Representation::__initialization(char *s) {
	if(!(the_global_discretization.u_h().matrix_representation()))
		the_global_discretization.u_h().matrix_representation() = this; // register mr for the variable
	// setup equation
	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] = new int[total_node_no*ndf];
	for(int i = 1; i < total_node_no; i++) eqn[i] = eqn[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_global_discretization.gh_on_gamma_h().gh_array()[i];
		if(gh(j) == gh_on_Gamma_h::Dirichlet) eqn[i][j] = Matrix_Representation::Fixed_Variable; // fixed degree of freedom
		else if(gh.tie_node_no() != -1) eqn[i][j] = Matrix_Representation::Tied_Node;            // tie node
		else eqn[i][j] = count++;
	}
   total_eqn_no = count;
   if(!s) { // default null string; use non-sparse matrix for rapid development of a small size problem
      the_lhs = new Global_Tensor(C0(total_eqn_no, total_eqn_no, (double*)0), *this);
      the_rhs = new Global_Tensor(C0(total_eqn_no, (double*)0), *this);
   } else { 
   	the_lhs = 0; the_rhs = 0; // initial state
	}
}

Matrix_Representation::Matrix_Representation(Global_Discretization& gd, char *s)
	: the_global_discretization(gd) {
   __initialization(s);
}

Matrix_Representation::Matrix_Representation(const Matrix_Representation& a) :
	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;
   if(a.the_lhs && the_lhs != a.the_lhs) {
   	delete the_lhs;
      int row_size = a.the_lhs->global_tensor().row_length(),
          col_size = a.the_lhs->global_tensor().col_length();
      the_lhs = new Global_Tensor(C0(row_size, col_size, (double*)0), *this);
      ((C0)*the_lhs) = ((C0)*(a.the_lhs));
   }
   if(((Matrix_Representation)a).rhs_ptr() && (((Matrix_Representation)(*this)).rhs_ptr() != ((Matrix_Representation)a).rhs_ptr()) ) {
   	delete &(rhs());
      int size = ((Matrix_Representation)a).rhs().global_tensor().length();
      ((Matrix_Representation)(*this)).rhs_ptr() = new Global_Tensor(C0(size, (double*)0), *this);
      ((C0)((Matrix_Representation)(*this)).rhs()) = ((C0)(((Matrix_Representation)a).rhs()));
   }
}

Matrix_Representation& Matrix_Representation::operator=(const Matrix_Representation& a) {
	if(this == &a) return *this;
	the_global_discretization = a.the_global_discretization;
   total_eqn_no = a.total_eqn_no;
   int total_node_no = a.the_global_discretization.omega_h().total_node_no(), // no distinct number
	    ndf = a.the_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;
   if(the_lhs) ((C0)*the_lhs) = ((C0)*(a.the_lhs));
   if(&(((Matrix_Representation)(*this)).rhs())) (C0)((Matrix_Representation)(*this)).rhs() = (C0)(((Matrix_Representation)a).rhs());
   return *this;
}

Matrix_Representation::~Matrix_Representation() {
	if(eqn[0] && eqn) {
   	delete [] eqn[0];
      eqn[0] = 0;
   }
   if(eqn) {
   	delete [] eqn;
      eqn = 0;
   }
   if(the_lhs) {
   	delete the_lhs;
      the_lhs = 0;
   }
   if(the_rhs) {
   	delete the_rhs;
      the_rhs = 0;
   }
}

Global_Tensor& Matrix_Representation::lhs() { return *the_lhs; }

Global_Tensor& Matrix_Representation::rhs() { return *the_rhs; }

U_h& U_h::operator=(C0& a) {
	if(!mr) {
		for(int i = 0; i < the_total_node_no; i++)
		for(int j = 0; j < ndf; j++)
			u_h_array()[i][j] = (double) a[i][j];
   } else {
		Global_Discretization &gd = mr->global_discretization();
		int** eqn = mr->equation_no_table();
		int ndf = gd.u_h()[0].no_of_dim();
		for(int i = 0; i < the_total_node_no; i++)
		for(int j = 0; j < ndf; j++) {
			if(eqn[i][j] >= 0) // free variable
				u_h_array()[i][j] = (double) a[eqn[i][j]];
			else if(eqn[i][j] == Matrix_Representation::Tied_Node) { // assign value of the tie node
      		int tie_node_no = u_h_array()[i].tie_node_no();
				u_h_array()[i][j] = (*this)[tie_node_no] [j];
      	}
      }
   }
   return *this;
}

U_h& U_h::operator+=(C0& a) {
	if(!mr) {
   	for(int i = 0; i < the_total_node_no; i++)
		for(int j = 0; j < ndf; j++)
			u_h_array()[i][j] += (double) a[i][j];
   } else {
		Global_Discretization gd = mr->global_discretization();
		int** eqn = mr->equation_no_table();
		int ndf = gd.u_h()[0].no_of_dim();
		for(int i = 0; i < the_total_node_no; i++)
		for(int j = 0; j < ndf; j++) {
			if(eqn[i][j] >= 0) // free variable
				u_h_array()[i][j] += (double) a[eqn[i][j]];
			else if(eqn[i][j] == Matrix_Representation::Tied_Node) { // a tie node can be free or fixed
      		int tie_node_no = u_h_array()[i].tie_node_no();
				u_h_array()[i][j] += (*this)[tie_node_no] [j];
      	}
      }
   }
   return *this;
}

U_h& U_h::operator-=(C0& a) {
	if(!mr) {
   	for(int i = 0; i < the_total_node_no; i++)
		for(int j = 0; j < ndf; j++)
			u_h_array()[i][j] -= (double) a[i][j];
   } else {
		Global_Discretization gd = mr->global_discretization();
		int** eqn = mr->equation_no_table();
		int ndf = gd.u_h()[0].no_of_dim();
		for(int i = 0; i < the_total_node_no; i++)
		for(int j = 0; j < ndf; j++) {
			if(eqn[i][j] >= 0) // free variable
				u_h_array()[i][j] -= (double) a[eqn[i][j]];
			else if(eqn[i][j] == Matrix_Representation::Tied_Node) { // a tie node can be free or fixed
      		int tie_node_no = u_h_array()[i].tie_node_no();
				u_h_array()[i][j] -= (*this)[tie_node_no] [j];
      	}
      }
   }
   return *this;
}

C0& Global_Tensor::operator+=(Element_Tensor& element_tensor) {
	int en = element_tensor.element_no(); 
   Omega_h& oh = element_tensor.matrix_representation().global_discretization().omega_h();
   Omega_eh elem = oh(en);
	int nen = elem.element_node_no(),
   	 first_node = elem[0],
	    ndf = element_tensor.matrix_representation().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 = element_tensor.matrix_representation().eqn[oh.node_order(node_no)][j];
            if(eq_no == Matrix_Representation::Tied_Node) {
      			int tie_node_no = element_tensor.matrix_representation().global_discretization().u_h().u_h_array()[oh.node_order(node_no)].tie_node_no();
               eq_no = element_tensor.matrix_representation().eqn[oh.node_order(tie_node_no)][j];
            }
            if(eq_no >=0 )
				the_global_tensor[eq_no] += element_tensor.element_tensor()[i*ndf+j];
         }
		}
	} else if(element_tensor.element_tensor().type() == matrix_type) {
		for(int i = 0; i < nen; i++)
		for(int j = 0; j < nen; j++) {
			int node_no1 = elem[i], node_no2 = elem[j];
			for(int k = 0; k < ndf; k++)
			for(int l = 0; l < ndf; l++) {
            int n1 = oh.node_order(node_no1), n2 = oh.node_order(node_no2),
                eqn1 = element_tensor.matrix_representation().eqn[n1][k],
 			    	 eqn2 = element_tensor.matrix_representation().eqn[n2][l];
            // converts tied node
            if(eqn1 == Matrix_Representation::Tied_Node) {
      			int tie_node_no = element_tensor.matrix_representation().global_discretization().u_h()[node_no1].tie_node_no();
               eqn1 = element_tensor.matrix_representation().eqn[oh.node_order(tie_node_no)][k];
            }
            if(eqn2 == Matrix_Representation::Tied_Node) {
      			int tie_node_no = element_tensor.matrix_representation().global_discretization().u_h()[node_no2].tie_node_no();
               eqn2 = element_tensor.matrix_representation().eqn[oh.node_order(tie_node_no)][l];
            }
            if(eqn1 >= 0 && eqn2 >= 0) // if none of the two dof is a fixed degree of freedom
					the_global_tensor[eqn1][eqn2] += element_tensor.element_tensor()[i*ndf+k][j*ndf+l];
			}