//=========================================================
// vs/fe output in FEMAP neutral file format
//=========================================================

void femap_material_property(ofstream& ofs1) {
    // neutral file header
    ofs1 << "   -1" << endl
         << "   100" << endl
         << "vs/fe output in FEMAP neutral file format" << endl
	 << "6.0," << endl
	 << "   -1" << endl;
    // layer data
    ofs1 << "   -1" << endl
         << "   413" << endl
	 << "1, 124," << endl
	 << "Default Layer" << endl
	 << "9999,124," << endl
	 << "Construction Layer" << endl
	 << "   -1" << endl;
    // materials
    ofs1 << "   -1" << endl
         << "   601" << endl
	 << "1,-601,104,0,0,1,0," << endl
	 << "Materials" << endl
	 << "10," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "25," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0," << endl
         << "200," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "50," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "60," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0,0,0," << endl
	 << "   -1" << endl;
    // properties
    ofs1 << "   -1" << endl
         << "   402" << endl
         << "1,110,1,17,1,0," << endl
         << "Properties" << endl
         << "0,0,0,0," << endl
         << "90," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0," << endl
         << "190," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
	 << "0.,0.,0.,0.,0.," << endl
         << "0," << endl
	 << "   -1" << endl;
}
void femap_node_and_element(Global_Discretization& gd, ofstream& ofs1) {
    Omega_h &oh = gd.omega_h();
    gh_on_Gamma_h &gh = gd.gh_on_gamma_h();
    // define node
    ofs1 << "   -1" << endl
         << "   403" << endl;
    for(int i = 0; i < oh.total_node_no(); i++)
    	ofs1 << oh[i].node_no()+1 // node ID starts from 1 in FEMAP
	     << ",0,0,1,46,"
	     << gh[i](0) << ","     // displ. constraints
	     << gh[i](1) 
	     << ",0,0,0,0,"                      // rot. constraints
	     << oh[i][0] << ","   // coordinates
	     << oh[i][1] << ",0.,0.," << endl;
    ofs1 << "   -1" << endl;
    // define element
    ofs1 << "   -1" << endl
         << "   404" << endl;
    for(i = 0; i < oh.total_element_no(); i++) {
    	ofs1 << oh(i).element_no()+1
	         << ",124,1,17,4,1,0,0,1,0,0,0," << endl;
	// element-node numbers
        for(int j = 0; j < oh(i).element_node_no(); j++)
           ofs1 << oh(i)[j]+1 << ",";
	    ofs1 << "0,0,0,0,0,0," << endl
	         << "0,0,0,0,0,0,0,0,0,0," << endl
	         << "0.,0.,0.," << endl
	         << "0.,0.,0.," << endl
	         << "0.,0.,0.," << endl
             << "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0," << endl;
    }
    ofs1 << "   -1" << endl;
}
void femap_constraint_and_load(Global_Discretization& gd, ofstream& ofs1) {
    Omega_h &oh = gd.omega_h();
    gh_on_Gamma_h &gh = gd.gh_on_gamma_h();
    // define constraint
    ofs1 << "   -1" << endl
         << "   506" << endl
		 << "1," << endl
	     << "Constraints" << endl;
    for(int i = 0; i < gh.gh_array().length(); i++) {
       ofs1 << gh.gh_array()[i].node_no()+1 // node ID starts from 1 in FEMAP
	        << ",120,1,"
	        << gh.gh_array()[i](0) << "," << gh.gh_array()[i](1) << ","
            << "0,0,0,0,0," << endl;
    }
    ofs1 << "-1,-1,-1,0,0,0,0,0,0,0," << endl
         << "-1,-1,-1,0,0,0,0,0,0,0," << endl
         << "-1,-1,-1,0,0,0,0,0,0,0," << endl
         << "-1,-1,-1,0,0,0,0,0,0,0," << endl
	     << "-1,-1,-1" << endl
	     << "   -1" << endl;
    // define load
    ofs1 << "   -1" << endl
         << "   507" << endl
		 << "1," << endl
	     << "Loadings" << endl
         << "0,0.,0,0,0," << endl
         << "0.,0.,0.," << endl
         << "0.,0.,0.," << endl
         << "0.,0.,0.," << endl
         << "0.,0.,0.," << endl
         << "0.,0.,0.,0," << endl
         << "0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.," << endl
         << "0,0,0," << endl
         << "0,0,0," << endl
         << "0.,0.,0.,0.," << endl
         << "0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.," << endl
         << "0.,0.,0.,0.,0.,0.,0.,0.,0.," << endl
         << "0,0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0," << endl
         << "0,0,0,0,0,0,0,0," << endl;
	for(i = 0; i < gh.gh_array().length(); i++) {
	   ofs1 << gh.gh_array()[i].node_no()+1; // node ID starts from 1 in FEMAP
	   if(gh.gh_array()[i](0) || gh.gh_array()[i](1))
	       ofs1 << ",3,10,1,0,0,0," << endl;
	   else ofs1 << ",1,10,1,0,0,0," << endl;
	   ofs1 << "1,1,1," << endl
	        << gh.gh_array()[i][0] << "," << gh.gh_array()[i][1] << ",0.0,0.0,0.0," << endl
            << "0,0,0,0,0," << endl
            << "0,0,0,0,0," << endl
            << "0,0,0," << endl
	        << "0.,0.,0.," << endl;
    }
    ofs1 << "-1,-1,-1,-1,-1,-1,0," << endl
         << "-1,-1,-1,-1,-1,-1,0," << endl
         << "-1,-1,-1,0.,0.,-1,1," << endl
         << "-1,-1,-1,0.,0.,-1,1," << endl
	 << "   -1" << endl;
}
void femap_view(Global_Discretization& gd, ofstream& ofs1) {
   Omega_h &oh = gd.omega_h();
   // define center for the view
   double min_x = 1.e20, min_y = 1.e20, max_x = -1.e20, max_y =-1.e20;
   for(int i = 0; i < oh.total_node_no(); i++) {
      if(oh[i][0] < min_x) min_x = oh[i][0];
      if(oh[i][0] > max_x) max_x = oh[i][0];
      if(oh[i][1] < min_y) min_y = oh[i][1];
      if(oh[i][1] > max_y) max_y = oh[i][1];
   }
   double center_x = (min_x+max_x)/2.0,
          center_y = (min_y+max_y)/2.0;
   ofs1 << "   -1" << endl
        << "   409" << endl
        << "1," << endl
        << "Default XY View" << endl
        << "0,0,0," << endl
        << "0.,0.,0.," << endl
        << center_x << "," << center_y << ",0.," << endl
	<< "1.,1.," << endl
        << center_x << "," << center_y << ",0.," << endl
        << "1.," << endl
        << "0.,0.,1.,1.," << endl
        << "2,0,1,1,0," << endl
        << "-1,-1,0,0,0,1,0,0,0,4000000," << endl
        << "9," << endl
        << "0,0,0," << endl
        << "0,0,0," << endl
        << "0,0,0," << endl
        << "0,0,0," << endl
        << "0,0,0," << endl
        << "0,0,0," << endl
        << "0,0,0," << endl
        << "0,0,0," << endl
        << "0,0,0," << endl
        << "100.,100.,1,7," << endl
        << "0,1,1.," << endl
        << "0.,0.,0.," << endl
        << "0.,0.,1.," << endl
        << "2.,1.,70.,0.5," << endl
        << "0.,25.,0.,0.,0.,0,100.,1000.,0.,0.,0.," << endl
        << "5.,90.,10.,10.,1.," << endl
        << "4,1,0,0,0,0,0,0,0,0,0.,0.,1.," << endl
        << "0,0,0,0,0,0,14,110," << endl
        << "0,1,1,1,0,1,0,1,1,0,1,1,0,1,1,1,0,0,1," << endl
        << "0,0,0.00000001,25.,100.,0.,0.,0.,20," << endl
        << "0,1,1,0,0,1,20.,0," << endl
        << "12," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,0.," << endl
        << "0,5,0,5,4,0.,25.," << endl
        << "4,16408,20,16504,100,16488," << endl
        << "0.,0.," << endl
        << "0.,0.,0.,0.," << endl
        << "9," << endl
        << "1.," << endl
        << "1.," << endl
        << "1.," << endl
        << "1.," << endl
        << "1.," << endl
        << "1.," << endl
        << "1.," << endl
        << "1.," << endl
        << "1.," << endl
        << "2," << endl
        << "<NULL>" << endl
        << "<NULL>" << endl
        << "0,0,0,0," << endl
        << "0.,0.,0.,0.," << endl
        << "0.,0.,0.,0.," << endl
        << "90,1,124,1,0," << endl
        << "0,60,0,1," << endl
        << "0,24,0,0," << endl
        << "0,120,0,0," << endl
        << "0,60,0,0," << endl
        << "0,24642,0,0," << endl
        << "0,124,0,0," << endl
        << "0,46,0,1," << endl
        << "0,120,0,1," << endl
        << "0,124,0,1," << endl
        << "0,124,0,0," << endl
        << "1,14,0,1," << endl
        << "0,62,0,0," << endl
        << "0,62,0,0," << endl
        << "0,10,0,1," << endl
        << "0,52,0,1," << endl
        << "0,4,0,1," << endl
        << "0,120,0,1," << endl
        << "0,12,0,1," << endl
        << "0,2,0,1," << endl
        << "0,120,0,1," << endl
        << "0,8312,0,1," << endl
        << "0,24600,0,0," << endl
        << "0,0,0,0," << endl
        << "1,123,0,1," << endl
        << "0,0,0,0," << endl
        << "2,124,0,1," << endl
        << "0,24636,0,0," << endl
        << "0,124,0,0," << endl
        << "0,4,0,0," << endl
        << "0,100,0,0," << endl
        << "0,124,0,1," << endl
        << "0,60,0,1," << endl
        << "0,56,0,1," << endl
        << "0,24,0,1," << endl
        << "0,8216,0,1," << endl
        << "0,4,0,0," << endl
        << "0,124,2,0," << endl
        << "0,0,1,1," << endl
        << "0,0,0,1," << endl
        << "1,124,5,1," << endl
        << "0,0,0,1," << endl
        << "0,24,0,1," << endl
        << "0,124,0,0," << endl
        << "0,100,0,1," << endl
        << "1,100,0,1," << endl
        << "0,0,0,1," << endl
        << "0,16,0,0," << endl
        << "0,124,4,1," << endl
        << "0,62,0,0," << endl
        << "2,124,1,1," << endl
        << "1,8254,0,0," << endl
        << "0,124,1,1," << endl
        << "1,0,5,1," << endl
        << "0,124,0,1," << endl
        << "0,100,0,1," << endl
        << "0,100,0,1," << endl
        << "1,46,0,1," << endl
        << "1,120,0,1," << endl
        << "1,4,0,1," << endl
        << "1,52,0,1," << endl
        << "1,24,0,1," << endl
        << "1,93,0,1," << endl
        << "1,12,0,1," << endl
        << "1,10,0,1," << endl
        << "1,104,0,1," << endl
        << "0,100,0,1," << endl
        << "0,24,0,1," << endl
        << "0,60,0,1," << endl
        << "0,104,0,1," << endl
        << "0,0,0,0," << endl
        << "0,90,0,1," << endl
        << "0,14,0,1," << endl
        << "0,8261,0,1," << endl
        << "0,19,0,0," << endl
        << "1,24,1,1," << endl
        << "0,4,0,0," << endl
        << "0,0,1,0," << endl
        << "0,0,0,0," << endl
        << "0,8201,0,1," << endl
        << "0,0,0,0," << endl
        << "0,0,1,1," << endl
        << "0,0,1,1," << endl
        << "0,0,1,1," << endl
        << "0,0,1,1," << endl
        << "0,0,1,1," << endl
        << "0,0,1,1," << endl
        << "0,0,1,1," << endl
        << "0,0,1,1," << endl
        << "0,0,1,1," << endl
        << "-1," << endl
        << "   -1" << endl
        << "   -1" << endl
        << "   412" << endl
        << "1,1,0,0," << endl
        << "   -1" << endl;
}
void femap_def(int time_step,Global_Discretization& gd, ofstream& ofs1) {
   Omega_h &oh = gd.omega_h();
   U_h &uh = gd.u_h();
	// define output set
	ofs1 << "   -1" << endl
		 << "   450" << endl
		 << time_step << endl
		 << "Step " << time_step << endl
		 << "0,1," << endl
		 << "0,0," << endl
		 << "1" << endl
		 << "<NULL>" << endl
		 << "   -1" << endl;
	// define deformation
	// 1. define total displacement; Vector ID = 1
	double *total_displacement, min = 1.e20, max = -1.e20;
	int min_id, max_id;
	total_displacement = new double[oh.total_node_no()];
	for(int i = 0; i < oh.total_node_no(); i++) {
		double ux = uh[i][0], uy = uh[i][1];
		total_displacement[i] = sqrt(ux*ux+uy*uy);
		if(total_displacement[i] < min) {
			min = total_displacement[i];
			min_id = i;
		}
		if(total_displacement[i] > max) {
			max = total_displacement[i];
			max_id = i;
		}
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",1,1," << endl
		 << "total displacement" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "2,3,4,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << total_displacement[i] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 2. define x-displacement; Vector ID = 2
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		if(uh[i][0] < min) { min = uh[i][0]; min_id = i; }
		if(uh[i][0] > max) { max = uh[i][0]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",2,1," << endl
		 << "x-displacement" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << uh[i][0] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 3. define y-displacement; Vector ID = 3
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		if(uh[i][1] < min) { min = uh[i][1]; min_id = i; }
		if(uh[i][1] > max) { max = uh[i][1]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",3,1," << endl
		 << "y-displacement" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << uh[i][1] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 4. define z-displacement; Vector ID = 4
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",4,1," << endl
		 << "z-displacement" << endl
		 << "0.0,0.0,0.0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "1,1,1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << ",0.," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
}
void femap_principal_stresses(int time_step, Global_Discretization& gd, ofstream& ofs1, C0& sig) {
	// compute principal stresses and principal directions
	double** principal_stress; // p1, p2, theta
	principal_stress = new double* [oh.total_node_no()];
	principal_stress[0] = new double [oh.total_node_no()*3];
	for(int i = 1; i < oh.total_node_no(); i++) 
		principal_stress[i] = principal_stress[i-1]+3;
	double mean, diff, rho;
	for(i = 0; i < oh.total_node_no(); i++) {
		mean = (sig[i][0]+sig[i][1])/2.0;
		diff = (sig[i][0]-sig[i][1])/2.0;
		rho = sqrt(diff*diff+(double)(sig[i][2]*sig[i][2]));
	    principal_stress[i][0] = mean + rho;
		principal_stress[i][1] = mean - rho;
		principal_stress[i][2] = 3.1416/4.0;
		if(fabs(diff) > 1.e-6) principal_stress[i][2] = atan(((double)sig[i][2])/diff)/2.0;
	}
	// define principal stresses and principal directions
	// 1. define maximum principal stresses; Vector ID = 5
	double min = 1.e20, max = -1.e20;
	int min_id, max_id;
	for(i = 0; i < oh.total_node_no(); i++) {
		if(principal_stress[i][0] < min) {
			min = principal_stress[i][0];
			min_id = i;
		}
		if(principal_stress[i][0] > max) {
			max = principal_stress[i][0];
			max_id = i;
		}
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",5,1," << endl
		 << "p1" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "6,7,8,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << principal_stress[i][0] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 2. define maximum principal stress x-component; Vector ID = 6
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		double x_component = principal_stress[i][0]*cos(principal_stress[i][2]);
		if(x_component < min) { min = x_component; min_id = i; }
		if(x_component > max) { max = x_component; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",6,1," << endl
		 << "p1: x-component" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << principal_stress[i][0]*cos(principal_stress[i][2]) << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 3. define y-displacement; Vector ID = 7
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		double y_component = principal_stress[i][0]*sin(principal_stress[i][2]);
		if(y_component < min) { min = y_component; min_id = i; }
		if(y_component > max) { max = y_component; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",7,1," << endl
		 << "p1: y-component" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << principal_stress[i][0]*sin(principal_stress[i][2]) << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 4. define z-displacement; Vector ID = 8
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",8,1," << endl
		 << "p1: z-component" << endl
		 << "0.0,0.0,0.0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "1,1,1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << ",0.," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 5. define minimun principal stresses; Vector ID = 11
	min = 1.e20, max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		if(principal_stress[i][1] < min) {
			min = principal_stress[i][1];
			min_id = i;
		}
		if(principal_stress[i][1] > max) {
			max = principal_stress[i][1];
			max_id = i;
		}
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",11,1," << endl
		 << "p2" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "12,13,14,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << principal_stress[i][1] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 6. define minimum principal stress x-component; Vector ID = 12
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		double x_component = principal_stress[i][1]*sin(principal_stress[i][2]);
		if(x_component < min) { min = x_component; min_id = i; }
		if(x_component > max) { max = x_component; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",12,1," << endl
		 << "p2: x-component" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << principal_stress[i][1]*sin(principal_stress[i][2]) << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 7. define minimun principal stress y-component; Vector ID = 13
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		double y_component = principal_stress[i][1]*cos(principal_stress[i][2]);
		if(y_component < min) { min = y_component; min_id = i; }
		if(y_component > max) { max = y_component; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",13,1," << endl
		 << "p2: y-component" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << principal_stress[i][1]*cos(principal_stress[i][2]) << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 8. define minimun principal stress z-component; Vector ID = 14
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",14,1," << endl
		 << "p2: z-component" << endl
		 << "0.0,0.0,0.0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "1,1,1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << ",0.," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
}
void femap_yield(int time_step, Global_Discretization& gd, ofstream& ofs1, C0& yield) {
	// define von Mises yield; Vector ID = 9
	double min = 1.e20, max = -1.e20;
	int min_id, max_id;
	for(int i = 0; i < oh.total_node_no(); i++) {
		if(((double)yield[i]) < min) { min = yield[i]; min_id = i; }
		if(((double)yield[i]) > max) { max = yield[i]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",9,1," << endl
		 << "von Mises yield" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << yield[i] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
}
void femap_sig(int time_step, Global_Discretization& gd, ofstream& ofs1, C0& sig) {
	// 1. define sig_xx; Vector ID = 25
	double min = 1.e20, max = -1.e20;
	int min_id, max_id;
	for(int i = 0; i < oh.total_node_no(); i++) {
		if(((double)sig[i][0]) < min) { min = sig[i][0]; min_id = i; }
		if(((double)sig[i][0]) > max) { max = sig[i][0]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",25,1," << endl
		 << "sig_xx" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << sig[i][0] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 2. define sig_yy; Vector ID = 26
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		if(((double)sig[i][1]) < min) { min = sig[i][1]; min_id = i; }
		if(((double)sig[i][1]) > max) { max = sig[i][1]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",26,1," << endl
		 << "sig_yy" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << sig[i][1] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 3. define sig_xy; Vector ID = 27
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		if(((double)sig[i][2]) < min) { min = sig[i][2]; min_id = i; }
		if(((double)sig[i][2]) > max) { max = sig[i][2]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",27,1," << endl
		 << "sig_xy" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << sig[i][2] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
}
void femap_eps(int time_step, Global_Discretization& gd, ofstream& ofs1, C0& eps) {
	// 1. define eps_xx; Vector ID = 28
	double min = 1.e20, max = -1.e20;
	int min_id, max_id;
	for(int i = 0; i < oh.total_node_no(); i++) {
		if(((double)eps[i][0]) < min) { min = eps[i][0]; min_id = i; }
		if(((double)eps[i][0]) > max) { max = eps[i][0]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",28,1," << endl
		 << "eps_xx" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << eps[i][0] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 2. define eps_yy; Vector ID = 29
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		if(((double)eps[i][1]) < min) { min = eps[i][1]; min_id = i; }
		if(((double)eps[i][1]) > max) { max = eps[i][1]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",29,1," << endl
		 << "eps_yy" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << eps[i][1] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 3. define eps_xy; Vector ID = 30
	min = 1.e20; max = -1.e20;
	for(i = 0; i < oh.total_node_no(); i++) {
		if(((double)eps[i][2]) < min) { min = eps[i][2]; min_id = i; }
		if(((double)eps[i][2]) > max) { max = eps[i][2]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << time_step << ",30,1," << endl
		 << "eps_xy" << endl
		 << min << "," << max << "," << fabs(max) << "," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << "0,0,0,0,0,0,0,0,0,0," << endl
		 << min_id+1 << "," << max_id+1 << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << i+1 << "," << eps[i][2] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
}
void femap_output_model(Global_Discretization& gd, ofstream& ofs1) {
   femap_material_property(ofs1);
   femap_node_and_element(gd, ofs1);
   femap_constraint_and_load(gd, ofs1);
   femap_view(gd, ofs1);
}
void femap_output_data(int time_step, Global_Discretization& gd, ofstream& ofs1, C0& sig, C0& eps, C0& yield) {
   // deformation plot and displacement component contours; Vector ID 1..4
   femap_def(time_step+1, gd, ofs1);
   // for stress and strain cross-hair plots; 
   // Vector ID 5..8 for maximun principal stress, Vector ID 11..14 for minimun principal stress
   // Vector ID 15..18 for maximun principal stress, Vector ID 21..24 for minimun principal stress
   femap_principal_stresses(time_step+1, gd, ofs1, sig);
   // component contours for von Mises yield; Vector ID 9
   femap_yield(time_step+1, gd, ofs1, yield);
   // component contours for stresses (Vector ID 25..27) and strains (Vector ID 28..30); 
   femap_sig(time_step+1, gd, ofs1, sig);
   femap_eps(time_step+1, gd, ofs1, eps);
}