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

Nodal_Constraint::Nodal_Constraint(int nn, int ndf, double* value, int* flag)
   : Nodal_Value(nn, ndf, value)  {
   constraint_flag = new int[ndf];
   if(flag) for(int i = 0; i < ndf; i++) constraint_flag[i] = flag[i];
   else for(int i = 0; i < ndf; i++) constraint_flag[i] = 0;
}

Nodal_Constraint::Nodal_Constraint(const Nodal_Constraint& a) : Nodal_Value(a) {
	int nd = no_of_dim();
	constraint_flag = new int[nd];
   for(int i = 0; i < nd; i++)
   	constraint_flag[i] = a.constraint_flag[i];
}

Nodal_Constraint& Nodal_Constraint::operator=(const Nodal_Constraint& a) {
	if(this != &a) {
   	Nodal_Value::operator=(a); // assign base class to this
		if(no_of_dim() != a.no_of_dim()) {
      	int nd = no_of_dim() = a.no_of_dim();
      	if(constraint_flag) delete [] constraint_flag;
         constraint_flag = new int[nd];
      }
      for(int i = 0; i < no_of_dim(); i++)
      	constraint_flag[i] = a.constraint_flag[i];
   }
   return (*this);
}

void gh_on_Gamma_h::__initialization(int df, Omega_h& oh) {
	total_node_no = oh.total_node_no();
   ndf = df;
   for(int i = 0; i < total_node_no; i++) {
      Node& node = oh.node_array()[i];
   	Nodal_Constraint *gh = new Nodal_Constraint(node.node_no(), ndf);
      gh->tie_node_no() = node.tie_node_no();
      the_gh_array.add(gh);
   }
}

Nodal_Constraint& gh_on_Gamma_h::operator[](int nn) {
	Dynamic_Array_Iterator<Nodal_Constraint> iter(the_gh_array);
   while(iter->node_no() != nn && iter.more()) iter++;
   return *(iter.current());
}

int gh_on_Gamma_h::node_order(int nn) {
	Dynamic_Array_Iterator<Nodal_Constraint> iter(the_gh_array);
   while(iter->node_no() != nn && iter.more()) iter++;
   return iter.index();
}

void U_h::__initialization(int df, Omega_h& oh) {
	ndf = df;
	mr = 0;
   the_total_node_no = oh.total_node_no();
   for(int i = 0; i < the_total_node_no; i++) {
      Node& node = oh.node_array()[i];
      Nodal_Value* uh = new Nodal_Value(node.node_no(), ndf);
      uh->tie_node_no() = node.tie_node_no();
      the_u_h_array.add(uh);
   }
}

Nodal_Value& U_h::operator[](int nn) {
	Dynamic_Array_Iterator<Nodal_Value> iter(the_u_h_array);
   while(iter->node_no() != nn && iter.more()) iter++;
   return *(iter.current());
}

int U_h::u_h_order(int nn) {
	Dynamic_Array_Iterator<Nodal_Value> iter(the_u_h_array);
   while(iter->node_no() != nn && iter.more()) iter++;
   return iter.index();
}

U_h& U_h::operator=(gh_on_Gamma_h& gh) {
   for(int i = 0; i < the_total_node_no; i++)
   for(int j = 0; j < ndf; j++)
      if(gh.gh_array()[i](j) == gh_on_Gamma_h::Dirichlet) {
      	int tnn = the_u_h_array[i].tie_node_no();
      	if(tnn == -1) the_u_h_array[i][j] = gh.gh_array()[i][j];
         //else (*this)[tnn][j] = gh.gh_array()[i][j];
         else the_u_h_array[i][j] = gh[tnn][j];
      }
   return (*this);
}

U_h::operator C0() {
	C0 ret_val;
	if(ndf == 1) { // return a vector
   	ret_val &= C0(the_total_node_no, (double*)0);
   	for(int i = 0; i < the_total_node_no; i++)
      		ret_val[i] = (*this)[i][0];
   } else { // return a matrix
   	ret_val &= C0(the_total_node_no, ndf, (double*)0);
   	for(int i = 0; i < the_total_node_no; i++)
   	for(int j = 0; j < ndf; j++)
      	ret_val[i][j] = the_u_h_array[i][j];
   }
   return ret_val;
}

ostream& operator <<(ostream& s, U_h& uh) {
	for(int i = 0; i < uh.total_node_no(); i++)
   	s << (uh.u_h_array()[i]) << endl;
   return s;
}

Global_Discretization& Global_Discretization::operator=(const Global_Discretization& a) {
	if(this != &a) {
   	the_omega_h = a.the_omega_h;
      the_u_h = a.the_u_h;
      the_gh_on_gamma_h = a.the_gh_on_gamma_h;
      type_id = a.type_id;
   }
   return *this;
}

C0 Global_Discretization::element_coordinate(int en) {
	int nen = the_omega_h(en).element_node_no(),
   	 nn  = the_omega_h(en).element_node_array()[0],
	    nsd = the_omega_h[nn].no_of_dim();
   C0 ret_val;
   if(nsd > 1) { // build a matrix
		ret_val &= C0(nen, nsd, (double*)NULL);
		for(int i = 0; i < nen; i++) {
			int node_no = the_omega_h(en)[i];
			Node& node = the_omega_h[node_no];
			for(int j = 0; j < nsd; j++)
				ret_val[i][j] = node[j];
		}
   } else { // build a vector
		ret_val &= C0(nen, (double*)NULL);
		for(int i = 0; i < nen; i++) {
			int node_no = the_omega_h(en)[i];
			Node& node = the_omega_h[node_no];
         ret_val[i] = node[0];
		}
   }
	return ret_val;
}

C0 Global_Discretization::element_fixed_variable(int en) {
	int nen = the_omega_h(en).element_node_no(),
   	 nn  = the_omega_h(en).element_node_array()[0],
	    ndf = the_u_h[nn].no_of_dim();
	C0 ret_val;
   ret_val &= C0(nen*ndf, (double*)NULL);
   for(int i = 0; i < nen; i++) {
   	int node_no = the_omega_h(en)[i];
      Nodal_Constraint& nc= the_gh_on_gamma_h[node_no];
      for(int j = 0; j < ndf; j++)
         //if(the_gh_on_gamma_h.gh_array()[node_no](j)
         if(the_gh_on_gamma_h[node_no](j)
            == gh_on_Gamma_h::Dirichlet) ret_val[i*ndf+j] = nc[j];
         else ret_val[i*ndf+j] = 0.0;
   }
	return ret_val;
}

C0 Global_Discretization::element_free_variable(int en) {
// residual = - (element_stiffness*element_free_variable)
	int nen = the_omega_h(en).element_node_no(),
   	 nn  = the_omega_h(en).element_node_array()[0],
	    ndf = the_u_h[nn].no_of_dim();
	C0 ret_val;
   ret_val &= C0(nen*ndf, (double*)NULL);
   for(int i = 0; i < nen; i++) {
      int node_no = the_omega_h(en)[i];;
      //Nodal_Value& var = the_u_h.u_h_array()[node_no];
      Nodal_Value& var = the_u_h[node_no];
      for(int j = 0; j < ndf; j++)
         //if(the_gh_on_gamma_h.gh_array()[node_no](j) == gh_on_Gamma_h::Neumann)
         if(the_gh_on_gamma_h[node_no](j) == gh_on_Gamma_h::Neumann)
            ret_val[i*ndf+j] = var[j];
   }
	return ret_val;
}

C0 Global_Discretization::element_nodal_force(int en) {
	int nen = the_omega_h(en).element_node_no(),
   	 nn  = the_omega_h(en).element_node_array()[0],
	    ndf = the_u_h[nn].no_of_dim();
   C0 ret_val;
   ret_val &= C0(nen*ndf, (double*)NULL);
   for(int i = 0; i < nen; i++) {
      int node_no = the_omega_h(en)[i];
      Nodal_Constraint& nc= the_gh_on_gamma_h[node_no];
      for(int j = 0; j < ndf; j++)
         //if(the_gh_on_gamma_h.gh_array()[node_no](j)
         if(the_gh_on_gamma_h[node_no](j)
            == gh_on_Gamma_h::Neumann) ret_val[i*ndf+j] = nc[j];
         else ret_val[i*ndf+j] = 0.0;
   }
	return ret_val;
}