omega_h.h

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

//==============================================================================
// 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 conceptual design class that does not really exist
//--USER PROGRAM----------------------------------------------------------------
//    Finite_Element_Approximation(gd); // Global_Discretization <--- Element_Formuoation
//                                      // eliminated by Automatic Type Register idiom
//------------------------------------------------------------------------------
//    reduce to top-down structure
//          Finite_Element_Approximation = { Global_Discretization, Element_Formulation}
//                       ^
//                       |
//            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
//
//==============================================================================

#ifndef __OMEGA_H_H
#define __OMEGA_H_H

//==============================================================================
// Part A. Definitions of Physical Entities--
//         classes: Node,
//                  Omega_eh (discretized element domain),
//                  Omega_h (discretized 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 Dynamic_Array<T> 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
//==============================================================================

// class Nodal_Value is an abstract class that Node and u_h derived from
class Nodal_Value {
protected:
	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);
   Nodal_Value(const Nodal_Value&);
   Nodal_Value& operator=(const Nodal_Value&);
	virtual ~Nodal_Value() { delete [] the_value; }
	int& node_no() { return the_node_no; }
   const int node_no() const { return the_node_no; }
	int& tie_node_no() { return the_tie_node_no; }
   const int tie_node_no() const { return the_tie_node_no; }
	int& no_of_dim() { return nd; }
   const int no_of_dim() const { 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 = 0, int nd = 0, double* coor = NULL) : Nodal_Value(nn, nd, coor) {}
   Node(const Node&);
   Node& operator=(const Node&);
	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 = 0, int etn = 0, int mtn = 0, int nen = 0, int* ena = NULL);
   Omega_eh(const Omega_eh&);
   virtual ~Omega_eh() { delete [] the_element_node_array; }
   Omega_eh& operator=(const Omega_eh&);
	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 { // DISCRETIZED GLOBAL DOMAIN
protected:
	Dynamic_Array<Node> the_node_array;
	Dynamic_Array<Omega_eh> the_omega_eh_array;
public:
	Omega_h(); // constructor left for customized problem definition
	Omega_h(int) { /* do nothing */ } // to be used as default value in Global_Discretization
   Omega_h(const Omega_h&);
   virtual ~Omega_h() {
   	the_node_array.Dynamic_Array<Node>::~Dynamic_Array<Node>();
      the_omega_eh_array.Dynamic_Array<Omega_eh>::~Dynamic_Array<Omega_eh>();
   }
   Dynamic_Array<Node>& node_array() { return the_node_array; }
	Node& operator[](int nn); // node selector
	Node& operator[](int nn) const; // node selector
	int node_order(int nn) const;
	int total_node_no() const { return the_node_array.length(); }
	int total_distinct_node_no() const; // total number of distinct nodes
   void set(Node*);

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