//=========================================================
// 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(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
		 << "1" << endl
		 << "Data Set 1 "<< 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
		 << "1,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 << "," << max_id << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << oh.node_array()[i].node_no() << "," << 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
		 << "1,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 << "," << max_id << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << oh.node_array()[i].node_no() << "," << 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
		 << "1,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 << "," << max_id << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << oh.node_array()[i].node_no() << "," << uh[i][1] << "," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
	// 4. define z-displacement; Vector ID = 4
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << "1,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 << oh.node_array()[i].node_no() << ",0.," << endl;
	ofs1 << "-1,0.," << endl
		 << "   -1" << endl;
}
void femap_J2(Global_Discretization& gd, ofstream& ofs1, C0& J2) {
    Omega_h &oh = gd.omega_h();
	// define J2; 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)J2[i]) < min) { min = J2[i]; min_id = i; }
		if(((double)J2[i]) > max) { max = J2[i]; max_id = i; }
	}
	ofs1 << "   -1" << endl
         << "   451" << endl
		 << "1,9,1," << endl
		 << "J2d" << 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 << "," << max_id << ",1,7," << endl
		 << "1,1,1," << endl;
	for(i = 0; i < oh.total_node_no(); i++) 
		 ofs1 << oh.node_array()[i].node_no() << "," << J2[i] << "," << 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(Global_Discretization& gd, ofstream& ofs1, C0& J2) {
   // deformation plot and displacement component contours; Vector ID 1..4
   femap_def(gd, ofs1);
   // component contours for von Mises yield; Vector ID 9
   femap_J2(gd, ofs1, J2);
}