matrix_representation.h

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

#ifndef __MATRIX_REPRESENTATION
#define __MATRIX_REPRESENTATION

//==============================================================================
// Part D. Matrix Representation
//       symbolic processing is hidden inside the solver of sparse matrix
//       Permutation matrices are used to pivoting and minimized profile
//       A x = b      =>       (P*A*Q)*Q^t*x = P*b,
//       P: row permutation matrix, Q: column permutation matrix
//       Q = P^t for symmetrical matrix to preserve symmetry
//==============================================================================
class Matrix_Representation {
	friend class Element_Tensor;
   friend class Global_Tensor;
   C0 the_global_nodal_value;
protected:
	Global_Tensor *the_lhs, *the_rhs;
	Global_Discretization &the_global_discretization;
	// eqn[node_no][ndf] stores equation number in sparse matrix,
	// therefore, eqn is a relational table
	int total_eqn_no, **eqn;
public:    
	void __initialization(char *s);
	enum assembly_switch { ALL = 1, // for LHS+RHS; default for member function assembly
                          LHS, RHS, REACTION,
                          STRESS, NODAL_STRESS, STRAIN, NODAL_STRAIN,
                          FLUX, NODAL_FLUX,      // for heat conduction and irrotational flow problems
                          GAUSS_SCALAR, NODAL_SCALAR}; // for yield ratio
   static assembly_switch Assembly_Switch;
   static int axisymmetric_flag;
	// profile/linked, sym/unsym sparse matrices
	Matrix_Representation(Global_Discretization& gd = Global_Discretization(), char* s = NULL);
	// set-up Sparse Matrix data structure,
	// set eqn[node_no][dof] = un-permutated equation number, or
	//                         = -1 if is a constrained dof
	//                         = -2 if is a tied node
   enum { Tied_Node = -2, Fixed_Variable };
   Matrix_Representation(const Matrix_Representation&);
   Matrix_Representation& operator=(const Matrix_Representation&);
	virtual ~Matrix_Representation();
   virtual Global_Tensor& lhs();
   virtual Global_Tensor& rhs();
   virtual Global_Tensor*& rhs_ptr() { return the_rhs; }
	virtual Global_Discretization& global_discretization() { return the_global_discretization; }
	virtual int** equation_no_table() { return eqn; }
	virtual void assembly(int nodal_load_flag = TRUE, int en1 = -1, int en2 = -1, int eni = 1);
	C0& global_nodal_value() { return the_global_nodal_value; }
};

class Global_Tensor {
	C0 the_global_tensor;
	Matrix_Representation &the_matrix_representation;
public:
	Global_Tensor(C0& gt = C0(), Matrix_Representation &mr = Matrix_Representation())
   	: the_global_tensor(gt), the_matrix_representation(mr) {}
   Global_Tensor(const Global_Tensor& a)
   	: the_matrix_representation(a.the_matrix_representation) {
   	the_global_tensor &= a.the_global_tensor; }
   virtual ~Global_Tensor() {}
   Global_Tensor& operator=(const Global_Tensor& a) {
   	if(this == &a) return *this;
   	the_matrix_representation = a.the_matrix_representation;
      the_global_tensor = a.the_global_tensor;
      return *this;
   }
   operator C0() { return the_global_tensor; }
   C0& global_tensor() { return the_global_tensor; }
   C0 global_tensor() const { return the_global_tensor; }
	C0& operator+=(Element_Tensor& element_tensor);
};

class Element_Tensor {
	int the_element_no;      // for lhs tensor is a sparse matrix, for rhs tensor is
	C0& the_element_tensor;
	Matrix_Representation &the_matrix_representation;
	friend class Global_Tensor;
public:
	Element_Tensor(int en, C0& element_tensor = C0(), Matrix_Representation &mr = Matrix_Representation())
   	: the_element_no(en), the_element_tensor(element_tensor), the_matrix_representation(mr) {}
   virtual ~Element_Tensor() {}
   operator C0() { return the_element_tensor; }
   Matrix_Representation& matrix_representation() { return the_matrix_representation; }
   int element_no() { return the_element_no; }
	C0& element_tensor() { return the_element_tensor; }
   C0 element_tensor() const { return the_element_tensor; }
};
#endif