global_discretization.h

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

#include <stdlib.h>
#include "include\vs.h"

//==============================================================================
// Part A. Definitions of Physical Entities--
//         classes: Node,
//                  Omega_eh (element domain),
//                  Omega_h (global domain)
//
// Choice of Data Structure of NODE, ELEMENT, U_H, ...etc.
// NODE can use SET or BAG data structures in Finite Element, elements in a
// SET is unique where elements in a BAG is not. It will be more flexible for
// Finite Element Input Specification if BAG is used to contain NODE elements;
// i.e., same physical node can be input with two different node numbers, where
// their coordinates are used to decide if the two nodes are the same or not
//==============================================================================

// a simplified bag data structure taken from Bruce Eckel "Black Belt C++" p.185-6
template<class T> class Bag_Iterator;
template<class T>
class Bag {
	T ** v;
   int size, next;
public:
	Bag(int ow = TRUE) : v(NULL), size(0), next(0) {}
   ~Bag();
   void add(T* P);
   const int length() { return next; }
   T& operator[](int i) { return *v[i]; }
   friend class Bag_Iterator<T>;
};

template<class T>
class Bag_Iterator {
	Bag<T>& bag;
   int the_index;
public:
	Bag_Iterator(Bag<T>& B) : bag(B), the_index(0) {}
   T* operator->() { return bag.v[the_index]; }
   T* top() { return bag.v[the_index = 0]; }
   T* end() { return bag.v[the_index = bag.length()-1]; }
   T* current() { return bag.v[the_index]; }
   int more() { if(the_index < bag.next-1) return TRUE; else return FALSE; }
   int index() { return the_index; }
   T* operator++(int) { return bag.v[the_index++]; }
   T* operator++() { the_index++; return bag.v[the_index]; }
   T* operator--(int) { T* tmp = bag.v[the_index]; the_index--; return tmp; }
   T* operator--() { return bag.v[--the_index]; }
   int length() { return bag.length(); }
};

// class Nodal_Value is an abstract class that Node and u_h derived from
class Nodal_Value {
	int the_node_no,
	    nd,         // number of dimension
	    the_tie_node_no;
	    // if two nodes have same coordinates the second node will be tied
	    // to the first, the first node number is stored in tie_node_no,
	    // if there is no tie node tie_node_no = -1
	double* the_value; // store the attribue of the node
public:
	Nodal_Value(int nn, int ndim, double* a = NULL);
	virtual ~Nodal_Value() { delete [] the_value; }
	int& node_no() { return the_node_no; }
	int& tie_node_no() { return the_tie_node_no; }
	int& no_of_dim() { return nd; }
	double* value() { return the_value; }
	double& operator[](int i) { return the_value[i]; }
   friend ostream& operator<<(ostream&, Nodal_Value&);
};

class Node : public Nodal_Value { // coordinates as attributes of the nodal value
public:                           
	Node(int nn, int nd, double* coor = NULL) : Nodal_Value(nn, nd, coor) {}
	int operator==(Node&);
	int operator!=(Node& a) { return !(operator==(a)); }
};

class Omega_eh { // Element Domain
	int the_element_no,
	    the_element_type_no,  // element type number
	    the_material_type_no, // material type number
	    the_element_node_no,  // number of node per element
       *the_element_node_array; // array of nodes can be non-unique nodes
public:
	Omega_eh(int en, int etn, int mtn, int nen, int* ena = NULL);
	int& element_no() { return the_element_no; }
	int& element_type_no() { return the_element_type_no; }
	int& material_type_no() { return the_material_type_no; }
	int& element_node_no() { return the_element_node_no; }
   int* element_node_array() { return the_element_node_array; }
	int& operator[](int i) { return the_element_node_array[i]; }
	int operator==(Omega_eh&);
	int operator!=(Omega_eh& a) {return !(operator==(a)); }
};

class Omega_h { // GLOBAL DOMAIN
	Bag<Node> the_node_array;
	Bag<Omega_eh> the_omega_eh_array;
public:
	Omega_h(); // constructor left for customized problem definition 
	Bag<Node>& node_array() { return the_node_array; }
	Node& operator[](int nn); // node selector
	int node_order(int nn);
	int total_node_no() { return the_node_array.length(); }
	int total_distinct_node_no(); // total number of distinct nodes

	Bag<Omega_eh>& omega_eh_array() { return the_omega_eh_array; }
	Omega_eh& operator()(int en); // element selector
	int element_order(int en);
	int total_element_no() { return the_omega_eh_array.length(); } // total number of distinct nodes
};

//==============================================================================
// Part B: classes of free and fixed variables
//         gh_on_Gamma_h (boundary conditions)
//         U_h            free variables
//==============================================================================
class Nodal_Constraint : public Nodal_Value {
	int* constraint_flag; // Neumann B.C. = 0; Dirichlet B.C. = 1
public:
	Nodal_Constraint(int nn, int ndf, double* value = NULL, int* flag = NULL);
	~Nodal_Constraint() { delete [] constraint_flag; }
	int& operator()(int i) { return constraint_flag[i]; }
};

class gh_on_Gamma_h { // Boundary Condition := { u_h = g and f = h on Gamma_h }
	int total_node_no;
	Bag<Nodal_Constraint> the_gh_array;
   void __initialization(Omega_h*);
public:
	enum { Neumann = 0, Dirichlet };
	// constructor for customized problem definition 
	gh_on_Gamma_h(Omega_h*);
	Bag<Nodal_Constraint>& gh_array() { return the_gh_array; }
	Nodal_Constraint& operator[](int node_no);
	int node_order(int nn);
};

//==============================================================================
// Object-Oriented Design and Analysis of Finite Element Method
//==============================================================================
// PHYSICAL ENTITIES include nodes and elements
// Node <---- Omega_eh; i.e., element depends on node
// Omega_h = { all Node, all Omega_eh }; global discretized domain is defined as
//    a set of all nodes and all elements.
// a finite element is defined as a set of {Omega_eh, u_eh, N_e} by P. Ciarlet, 1976
//
// Dependency Graph is a complicated network
//                             (seperator)
//                    Omega_h  <-----> Element_Formulation
//                 ^         ^            ^          |
//               /            \          /           |
//             U_h <-----> Matrix_Representation     |
//              ^ (seperator)                        |
//              +------------------------------------+
// 
// forward declarations for linearized seperate compilation among Omega_h,
// U_h(u_h; g_h, h_h), Element_Formulation, and Matrix_Representation
// classes which form an unfortunate network,
// CYCLIC class dependancy GRAPH-- clique shows up
//
//  First select and preserve the relation: Omega_h <---- u_h,
//  that is the declaration of u_h will have to refer to Omega_h as
//  Omega_h oh;
//  U_h uh(oh);
//
// A. Global_Discretization = { Omega_h, U_h, gh_on_Gamma_h} is a
//    strong component design class
//            Global_Discretization = { Omega_h, u_h}
//                         ^
//                         |
//                 Element_Formulation
//    reduce three elements relation to to-down structure by using an int number:
//    melement_type_no to refer to Element_Formulation, and use Automatic Type
//    Register idiom to register more than one Element_Formulation
//--USER PROGRAM----------------------------------------------------------------
//    Global_Discretization gd(Omega_h, u_h, gh_on_Gamma_h)
//    Element_Formulation(int element_no, Global_Discretization& gd);
//    Element_Formulation* Element_Formulation::type_list = 0;
//    static HeatQ4 heatq4_instance(Element_Type_Register);
//------------------------------------------------------------------------------
//==============================================================================
// B. Finite_Element_Approximation = { Omega_h, u_h, Element_Formulation}
//                                 = { Global_Discretization, Element_Formulation}
//    is a strong component design class
//--USER PROGRAM----------------------------------------------------------------
//    Finite_Element_Approximation(gd); // Global_Discretization <--- Element_Formuoation
//                                      // eliminated by Automatic Type Register idiom
//------------------------------------------------------------------------------
//    reduce to top-down structure
//          Finite_Element_Approximation
//                       ^
//                       |
//            Matrix_Representation
//    the compliation dependancy is seperated between U_h and Matrix_Representation
//    cyclic dependancy by first declare a pointer to Matrix_Representation in U_h,
//    then use a a private function: Matrix_Representation::__initialization()
//    called by the constructor of the Matrix_Representation, that sets up the
//    pointer to Matrix_Representation in U_h to refer to the current instance of the
//    Matrix_Representation
//
//==============================================================================
class Matrix_Representation;
class U_h { // Nodal Variable is the "Basis" of finite element approximation
	int the_total_node_no, ndf;
	Matrix_Representation *mr;
	Bag<Nodal_Value> the_u_h_array;
   void __initialization(Omega_h*);
public:
	U_h(Omega_h* oh) { __initialization(oh); }
	Bag<Nodal_Value>& u_h_array() { return the_u_h_array; }
	Nodal_Value& operator[](int nn);
	int u_h_order(int nn);
   int total_node_no() { return the_total_node_no; }
	U_h& operator=(gh_on_Gamma_h& gh);
	// to be used in constructor of Matrix_Representation to close cyclic dependency
	Matrix_Representation* matrix_representation() { return mr; }
	U_h& operator=(C0&);
	U_h& operator+=(C0&);
	U_h& operator-=(C0&);
   friend ostream& operator<<(ostream&, U_h&);
};

//==============================================================================
//  an integrated class of Part A and B
//  --a "strong component design class" in Graph Theory
//==============================================================================
class Global_Discretization {
	Omega_h *the_omega_h;
	U_h *the_u_h;
	gh_on_Gamma_h *the_gh_on_gamma_h;
public:
	Global_Discretization(Omega_h* oh, gh_on_Gamma_h* gh, U_h* uh) :
		the_omega_h(oh), the_gh_on_gamma_h(gh), the_u_h(uh) {}
	Omega_h* omega_h() { return the_omega_h; }
	U_h& u_h() { return *the_u_h; }
	gh_on_Gamma_h& gh_on_gamma_h() { return *the_gh_on_gamma_h; }
	C0 element_coordinate(int en);
   // integrated-information services to Element_Formulation
	C0 element_fixed_variable(int en);
	C0 element_free_variable(int en);
   C0 element_nodal_force(int en);
   int operator==(const Global_Discretization&) const;
};