/**********************************************************************************
 * Copyright 2010 Christoph Pflaum 
 * 		Department Informatik Lehrstuhl 10 - Systemsimulation
 *		Friedrich-Alexander Universität Erlangen-Nürnberg
 * 
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 **********************************************************************************/

#include "../mympi.h"
#include "../abbrevi.h"
#include "../parameter.h"
#include "../math_lib/math_lib.h"
#include "../basics/basic.h"
#include "elements.h"
#include "../grid/parti.h"
#include "ug.h"
#include "marker.h"

#include <algorithm>
#include <iomanip>

//////////////////////////////////////////////////////////////////
// Object_based_on_ug member functions
//////////////////////////////////////////////////////////////////

Object_based_on_ug::Object_based_on_ug(Unstructured_grid* ug) {
  ug->Add_object(this);
}

void Object_based_on_ug::Unregister(Unstructured_grid* ug) {
  ug->Remove_object(this);
}


//////////////////////////////////////////////////////////////////
// other member functions
//////////////////////////////////////////////////////////////////

void Unstructured_grid::Add_object(Object_based_on_ug *obj) { 
#if 0
  std::ios_base::fmtflags f = cout.flags();
  cout << "Unstructured_grid::Add_object() obj = 0x" 
    << std::hex
    << obj << endl;
  cout.setf(f);
  //*static_cast<int*>(0)=1;
#endif
    all_objects_based_on_ug.push_back(obj);
}

void Unstructured_grid::Remove_object(Object_based_on_ug *obj) {
#if 0
  std::ios_base::fmtflags f = cout.flags();
  cout << "Unstructured_grid::Remove_object() obj = 0x" 
    << std::hex 
    << obj << endl;
  cout.setf(f);
#endif

  Object_list::iterator it;
  it=find(all_objects_based_on_ug.begin(),
    all_objects_based_on_ug.end(),
    obj);
  if(it != all_objects_based_on_ug.end())
    all_objects_based_on_ug.erase(it);
}

Not_constant_direction_marker::Not_constant_direction_marker(Unstructured_grid* grid) {
  num_hexahedra = grid->Give_number_hexahedra();

  not_const_x = new bool[num_hexahedra];
  not_const_y = new bool[num_hexahedra];
  not_const_z = new bool[num_hexahedra];

  for(int i=0;i<num_hexahedra;++i) {
    not_const_x[i] = true;
    not_const_y[i] = true;
    not_const_z[i] = true;
  }
}

Not_constant_direction_marker::~Not_constant_direction_marker()
{
    delete[] not_const_x;
    delete[] not_const_y;
    delete[] not_const_z;
}

void Not_constant_direction_marker::Mark_x_as_not_constant_at(int i, bool bb) {
  if(i<0 || i>=num_hexahedra) 
    cout << " error in Constant_direction_marker::Mark_x_as_const_at" << endl;
  else not_const_x[i] = bb;
}

void Not_constant_direction_marker::Mark_y_as_not_constant_at(int i, bool bb) {
  if(i<0 || i>=num_hexahedra) 
    cout << " error in Constant_direction_marker::Mark_y_as_const_at" << endl;
  else not_const_y[i] = bb;
}

void Not_constant_direction_marker::Mark_z_as_not_constant_at(int i, bool bb) {
  if(i<0 || i>=num_hexahedra) 
    cout << " error in Constant_direction_marker::Mark_z_as_const_at" << endl;
  else not_const_z[i] = bb;
}



D3vector transform_nothing(double x, double* gd) {
  return D3vector(0.0,0.0,0.0);
}

D3vector transform_nothing(double x, double y) {
  return D3vector(0.0,0.0,0.0);
}

Unstructured_grid::Unstructured_grid() {
#if 0
  cout << "Unstructured_grid::Unstructured_grid()" << endl;
#endif
  
  constructionParameters.resize(0);

  construction_hexahedron_points_done_yn = false;
  construction_done_yn                   = false;
  periodic                               = false;

  pointer_global_data = NULL;
  size_pointer_global_data = 0;
   
  num_points = 0;
  num_edges = 0;
  num_quadrangles = 0;
  num_hexahedra  = 0;

  points = NULL;
  edges  = NULL;
  quadrangles = NULL;
  hexahedra   = NULL;

  number_degree_of_freedom = 0;

  all_points = NULL;
  
  typFuerSlice = noWithZTyp;
}

Unstructured_grid::~Unstructured_grid() {
  Object_list::iterator iter;
  for(iter = all_objects_based_on_ug.begin();
      iter !=  all_objects_based_on_ug.end();
      ++iter)
  {
#if 0
  std::ios_base::fmtflags f = cout.flags();
  cout << "Unstructured_grid::~Unstructured_grid() Delete_data() obj = 0x" 
    <<  std::hex
    << (*iter) << endl;
  cout.setf(f);
#endif
    (*iter)->Delete_data();
  }
  all_objects_based_on_ug.clear();

  if(points != NULL)
      delete [] points;
  points = NULL;
  if(edges != NULL)
      delete [] edges;
  edges = NULL;
  if(quadrangles != NULL)
      delete [] quadrangles;
  quadrangles = NULL;
  if(hexahedra != NULL)
      delete [] hexahedra;
  hexahedra = NULL;

  delete [] pointer_global_data;

  delete all_points;

  delete zordering;
  if (my_not_constant_directions != NULL)
    delete my_not_constant_directions;
}

void Unstructured_grid::Set_number_points(int num_points_) {
  num_points = num_points_;
  if(developer_version)
     if(points!=NULL)
       cout << " error in Unstructured_grid::Set_number_points !" << endl;
  points = new Point_el[num_points];
}

void Unstructured_grid::Set_coordinate_point(int id, D3vector coord) {
  if(developer_version)
     if(points==NULL)
        cout << " error in Unstructured_grid::Set_coordinate_point " << endl; 
  points[id].Set_coordinate(coord);
}

void Unstructured_grid::Set_number_hexahedra(int num_hexahedra_) {
  // Set informations of hexahedra
  num_hexahedra = num_hexahedra_;
  if(developer_version)
     if(hexahedra!=NULL) 
        cout << " error in Unstructured_grid::Set_number_hexahedra !" << endl;
  hexahedra = new Hexahedron_el[num_hexahedra];
}

void Unstructured_grid::Set_hexahedron(int id,
				       int i_WSD, int i_ESD,
				       int i_WND, int i_END,
				       int i_WST, int i_EST,
				       int i_WNT, int i_ENT) {
  if(developer_version)
     if(construction_hexahedron_points_done_yn==true)
        cout << " error in Unstructured_grid::Set_hexahedron !" << endl;

  // Anwender muss sehr schlau sein beim Anwenden dieser
  // Funktion, da sonst Orientierung nicht stimmt.
//  std::cout <<
//   i_WSD << " " <<  i_ESD
//  << " " <<  i_WND << " " <<  i_END
//  << " " <<  i_WST << " " <<  i_EST
//  << " " <<  i_WNT << " " <<  i_ENT << std::endl;
  hexahedra[id].Set_id_corner(WSDdir3D,i_WSD);
  hexahedra[id].Set_id_corner(ESDdir3D,i_ESD);
  hexahedra[id].Set_id_corner(WNDdir3D,i_WND);
  hexahedra[id].Set_id_corner(ENDdir3D,i_END);
  hexahedra[id].Set_id_corner(WSTdir3D,i_WST);
  hexahedra[id].Set_id_corner(ESTdir3D,i_EST);
  hexahedra[id].Set_id_corner(WNTdir3D,i_WNT);
  hexahedra[id].Set_id_corner(ENTdir3D,i_ENT);


  
  hexahedra[id].Set_coord(WSDdir3D,points[i_WSD].Give_coordinate());
  hexahedra[id].Set_coord(ESDdir3D,points[i_ESD].Give_coordinate());
  hexahedra[id].Set_coord(WNDdir3D,points[i_WND].Give_coordinate());
  hexahedra[id].Set_coord(ENDdir3D,points[i_END].Give_coordinate());
  hexahedra[id].Set_coord(WSTdir3D,points[i_WST].Give_coordinate());
  hexahedra[id].Set_coord(ESTdir3D,points[i_EST].Give_coordinate());
  hexahedra[id].Set_coord(WNTdir3D,points[i_WNT].Give_coordinate());
  hexahedra[id].Set_coord(ENTdir3D,points[i_ENT].Give_coordinate());
}

void Unstructured_grid::Set_hexahedron(int id,
				       int i_WSD, int i_ESD,
				       int i_WND, int i_END,
				       int i_WST, int i_EST,
				       int i_WNT, int i_ENT,
				       D3vector coord_WSD, D3vector coord_ESD,
				       D3vector coord_WND, D3vector coord_END,
				       D3vector coord_WST, D3vector coord_EST,
				       D3vector coord_WNT, D3vector coord_ENT) {
  if(developer_version)
     if(construction_hexahedron_points_done_yn==true)
        cout << " error in Unstructured_grid::Set_hexahedron !" << endl;


  hexahedra[id].Set_id_corner(WSDdir3D,i_WSD);
  hexahedra[id].Set_id_corner(ESDdir3D,i_ESD);
  hexahedra[id].Set_id_corner(WNDdir3D,i_WND);
  hexahedra[id].Set_id_corner(ENDdir3D,i_END);
  hexahedra[id].Set_id_corner(WSTdir3D,i_WST);
  hexahedra[id].Set_id_corner(ESTdir3D,i_EST);
  hexahedra[id].Set_id_corner(WNTdir3D,i_WNT);
  hexahedra[id].Set_id_corner(ENTdir3D,i_ENT);


  hexahedra[id].Set_coord(WSDdir3D,coord_WSD);
  hexahedra[id].Set_coord(ESDdir3D,coord_ESD);
  hexahedra[id].Set_coord(WNDdir3D,coord_WND);
  hexahedra[id].Set_coord(ENDdir3D,coord_END);
  hexahedra[id].Set_coord(WSTdir3D,coord_WST);
  hexahedra[id].Set_coord(ESTdir3D,coord_EST);
  hexahedra[id].Set_coord(WNTdir3D,coord_WNT);
  hexahedra[id].Set_coord(ENTdir3D,coord_ENT);
}

elementTyp Give_corresponding_el_typ(elementTyp el_typ){
  if(el_typ == hexahedronEl)
     return quadrangleEl;
  if(el_typ == quadrangleEl)
     return edgeEl;
  if(developer_version)
     if(el_typ != edgeEl) cout << " errror in Give_corresponding_el_typ!" << endl;
  return pointEl;
}

void Calculate_corresponding_corners(dir3D_sons corner, dir3D_sons& SED_corner,
                                     dir3D_sons& NWD_corner, dir3D_sons& SWT_corner){
  if(corner == WSDdir3D){
     SED_corner = ESDdir3D;
     NWD_corner = WNDdir3D;
     SWT_corner = WSTdir3D;
  }
  else if(corner == ESDdir3D){
     SED_corner = ENDdir3D;
     NWD_corner = WSDdir3D;
     SWT_corner = ESTdir3D;
  }
  else if(corner == WNDdir3D){
     SED_corner = WSDdir3D;
     NWD_corner = ENDdir3D;
     SWT_corner = WNTdir3D;
  }
  else if(corner == ENDdir3D){
     SED_corner = WNDdir3D;
     NWD_corner = ESDdir3D;
     SWT_corner = ENTdir3D;
  }
  else if(corner == WSTdir3D){
     SED_corner = WNTdir3D;
     NWD_corner = ESTdir3D;
     SWT_corner = WSDdir3D;
  }
  else if(corner == ESTdir3D){
     SED_corner = WSTdir3D;
     NWD_corner = ENTdir3D;
     SWT_corner = ESDdir3D;
  }
  else if(corner == WNTdir3D){
     SED_corner = ENTdir3D;
     NWD_corner = WSTdir3D;
     SWT_corner = WNDdir3D;
  }
  else{
     if(developer_version)
        if(corner != ENTdir3D) cout << "error in Calculate_corresponding_corners!" << endl;
     SED_corner = ESTdir3D;
     NWD_corner = WNTdir3D;
     SWT_corner = ENDdir3D;
  }
}

void Calculate_corresponding_edges(dir3D_sons corner, Edges_cell& SD_edge,
                                   Edges_cell& WD_edge, Edges_cell& SW_edge){
  if(corner == WSDdir3D){
     SD_edge = SDed;
     WD_edge = WDed;
     SW_edge = SWed;
  }
  else if(corner == ESDdir3D){
     SD_edge = EDed;
     WD_edge = SDed;
     SW_edge = SEed;
  }
  else if(corner == WNDdir3D){
     SD_edge = WDed;
     WD_edge = NDed;
     SW_edge = NWed;
  }
  else if(corner == ENDdir3D){
     SD_edge = NDed;
     WD_edge = EDed;
     SW_edge = NEed;
  }
  else if(corner == WSTdir3D){
     SD_edge = WTed;
     WD_edge = STed;
     SW_edge = SWed;
  }
  else if(corner == ESTdir3D){
     SD_edge = STed;
     WD_edge = ETed;
     SW_edge = SEed;
  }
  else if(corner == WNTdir3D){
     SD_edge = NTed;
     WD_edge = WTed;
     SW_edge = NWed;
  }
  else{
     if(developer_version)
        if(corner != ENTdir3D) cout << "error in Calculate_corresponding_edges!" << endl;
     SD_edge = ETed;
     WD_edge = NTed;
     SW_edge = NEed;
  }
}

void Unstructured_grid::construction_hexahedron_points_done() {
  int num, num_q, ec, fc, c, i, e, q, ed, ismallest, smallest, number_neighbors, number_neighbors_hex, num_hex;
  int idA, idB, id, id_ed, id_quad, id_edge, id_point, id_corner;
  int id4[4];
  int idt[4];
  sorted_quad SQuad;
  sorted_pair SPair;
  unsorted_pair USPair;
  int c_SW, c_SE, c_NW, c_NE;
  int ed_W, ed_E, ed_S, ed_N;
  int id_search_quad, id_search_edge;
  int W_corner, E_corner;
  elementTyp neighbor_in_typ;
  int neighbor_in_id, id_global_ND, id_global_ST, id_global_SD, id_global_WD, id_global_SW,
      id_global_S, id_global_D, id_global_W;
  Edges_cell ed_cell, SD_edge, WD_edge, SW_edge;
  dir3D_sons corner, SED_corner, NWD_corner, SWT_corner;


  construction_hexahedron_points_done_yn = true;

  // allocate my_constant_directions
  my_not_constant_directions = new Not_constant_direction_marker(this);


  // Set informations of hexahedra, count edges and quadrangles and put in Topf
  for(num=0;num<num_hexahedra;++num) {
    // iterate edges
    for(ec=0;ec<12;++ec) {
      idA = hexahedra[num].Give_id_corner(Transform((Edges_cell)ec,Ldir1D));
      idB = hexahedra[num].Give_id_corner(Transform((Edges_cell)ec,Rdir1D));

      if(idB<idA)
      {
          id=idB;   idB=idA;   idA=id; hexahedra[num].orientation[ec] = false;
      }
      else
      {
          hexahedra[num].orientation[ec] = true;
      }

      topf_edges.put(idA,idB);
      hexahedra[num].Set_id_edge((Edges_cell)ec, topf_edges.Give_id(idA,idB));
    }
    // iterate quadrangles
    for(fc=0;fc<6;++fc) {
      for(i=0;i<4;++i) {
	id4[i] = hexahedra[num].Give_id_corner(corner_sur((dir3D)fc,
							  (dir2D_sons)i));
      }
      // find smallest corner id
      smallest  = id4[0];
      ismallest = 0;
      for(i=1;i<4;++i) {
        idt[i]=0; //mit 0 initialisieren
	if(id4[i] < smallest) { smallest = id4[i]; ismallest = i; }
      }

      if(SWdir2D==ismallest) {
	idt[SWdir2D] = id4[SWdir2D];
	idt[NEdir2D] = id4[NEdir2D];
	if(id4[NWdir2D] < id4[SEdir2D]) {
	  idt[NWdir2D] = id4[SEdir2D];
	  idt[SEdir2D] = id4[NWdir2D];
	}
	else {
	  idt[NWdir2D] = id4[NWdir2D];
	  idt[SEdir2D] = id4[SEdir2D];
	}
      }

      if(SEdir2D==ismallest) {
	idt[SWdir2D] = id4[SEdir2D];
	idt[NEdir2D] = id4[NWdir2D];
	if(id4[SWdir2D] < id4[NEdir2D]) {
	  idt[NWdir2D] = id4[NEdir2D];
	  idt[SEdir2D] = id4[SWdir2D];
	}
	else {
	  idt[NWdir2D] = id4[SWdir2D];
	  idt[SEdir2D] = id4[NEdir2D];
	}
      }

      if(NWdir2D==ismallest) {
	idt[SWdir2D] = id4[NWdir2D];
	idt[NEdir2D] = id4[SEdir2D];
	if(id4[NEdir2D] < id4[SWdir2D]) {
	  idt[NWdir2D] = id4[SWdir2D];
	  idt[SEdir2D] = id4[NEdir2D];
	}
	else {
	  idt[NWdir2D] = id4[NEdir2D];
	  idt[SEdir2D] = id4[SWdir2D];
	}
      }

      if(NEdir2D==ismallest) {
	idt[SWdir2D] = id4[NEdir2D];
	idt[NEdir2D] = id4[SWdir2D];
	if(id4[NWdir2D] < id4[SEdir2D]) {
	  idt[NWdir2D] = id4[SEdir2D];
	  idt[SEdir2D] = id4[NWdir2D];
	}
	else {
	  idt[NWdir2D] = id4[NWdir2D];
	  idt[SEdir2D] = id4[SEdir2D];
	}
      }
      topf_quadrangles.put(idt[SWdir2D],idt[SEdir2D],idt[NWdir2D],idt[NEdir2D]);

      hexahedra[num].Set_id_quadrangle((dir3D)fc,
		     topf_quadrangles.Give_id(idt[SWdir2D],idt[SEdir2D],
					      idt[NWdir2D],idt[NEdir2D]));
      }
  }

  // construct quadrangles and set informations of quadrangles
  num_quadrangles = topf_quadrangles.size();
  quadrangles = new Quadrangle_el[num_quadrangles];

  for(num=0;num<num_quadrangles;++num) {
    SQuad = topf_quadrangles.find(num);

    quadrangles[num].Set_id_corner(SWdir2D,SQuad.first);
    quadrangles[num].Set_id_corner(SEdir2D,SQuad.second);
    quadrangles[num].Set_id_corner(NWdir2D,SQuad.third);
    quadrangles[num].Set_id_corner(NEdir2D,SQuad.forth);

    
    // set edges
    for(ed=0;ed<4;++ed) {
      idA = quadrangles[num].Give_id_corner(Transform((dir2D)ed,Ldir1D));
      idB = quadrangles[num].Give_id_corner(Transform((dir2D)ed,Rdir1D));

      if(idB<idA)
      {
          id=idB;   idB=idA;   idA=id;
      }


      quadrangles[num].Set_id_edge((dir2D)ed, topf_edges.Give_id(idA,idB));
    }
  }

 
  
  // construct edges and set informations of edges
  num_edges = topf_edges.size();
  edges = new Edge_el[num_edges];

  for(num=0;num<num_edges;++num) {
    SPair = topf_edges.find(num);
    idA = SPair.first;
    idB = SPair.second;

    if(idB<idA) { id=idB;   idB=idA;   idA=id; }
    edges[num].Set_id_corners(idA,idB);
  }

  // construct neighbors of quadrangles
  for(num=0;num<num_hexahedra;++num){
      for(q=0; q<6; ++q){
          id_quad = hexahedra[num].Give_id_quadrangle((dir3D)q);
//           if ((dir3D)q == Tdir3D) cout << "top quad with num : " << id_quad << endl;
//         if ((dir3D)q == Ddir3D) cout << "Down quad with num : " << id_quad << endl;
          if(quadrangles[id_quad].Give_exists_exterior() == true){
             // Exists exterior neighbor false (default = true)
             quadrangles[id_quad].Set_exists_exterior(false);
             // Set interior volume
             quadrangles[id_quad].Set_id_interior_vol(num);
             quadrangles[id_quad].Set_typ_interior_vol(hexahedronEl);

             // Set interior faces
             // a) id of edges of Quad
             ed_W = quadrangles[id_quad].Give_id_edge(Wdir2D);
             ed_E = quadrangles[id_quad].Give_id_edge(Edir2D);
             ed_S = quadrangles[id_quad].Give_id_edge(Sdir2D);
             ed_N = quadrangles[id_quad].Give_id_edge(Ndir2D);
             // b) Iterate over all faces of hexahedron
             for(i=0;i<6;i++){
                 id_search_quad = hexahedra[num].Give_id_quadrangle((dir3D)i);
                 if(id_search_quad != id_quad){  // face is not the considered quad
                    for(ed = 0;ed<4;++ed){       // Iterate over all edges of faces
                                                  // and compare with edges of quad
                       id_edge = quadrangles[id_search_quad].Give_id_edge((dir2D)ed);
                       if(id_edge == ed_W)
                           quadrangles[id_quad].Set_id_interior_face(Wdir2D,id_search_quad);
                       else if(id_edge == ed_E)
                           quadrangles[id_quad].Set_id_interior_face(Edir2D,id_search_quad);
                       else if(id_edge == ed_S)
                           quadrangles[id_quad].Set_id_interior_face(Sdir2D,id_search_quad);
                       else if(id_edge == ed_N)
                           quadrangles[id_quad].Set_id_interior_face(Ndir2D,id_search_quad);
                     }
                  }
             }
             // c) Set type of interior faces
             for(i = 0;i<4;++i){
                 quadrangles[id_quad].Set_typ_interior_face((dir2D)i, quadrangleEl);
             }
             // Set interior edges
             // a) id of corners of Quad
             c_SW = quadrangles[id_quad].Give_id_corner(SWdir2D);
             c_SE = quadrangles[id_quad].Give_id_corner(SEdir2D);
             c_NW= quadrangles[id_quad].Give_id_corner(NWdir2D);
             c_NE = quadrangles[id_quad].Give_id_corner(NEdir2D);
             // b) Iterate over all edges of hexahedron
             for(i=0;i<12;i++){
                 id_search_edge = hexahedra[num].Give_id_edge((Edges_cell)i);
                 if(id_search_edge != ed_W &&
                     id_search_edge != ed_E &&
                     id_search_edge != ed_S &&
                     id_search_edge != ed_N){  // edge belongs not to the considered quad
                     for(c = 0;c<2;++c){      // Iterate over all corners of edges
                                               // and compare with corners of quad
                       id_point = edges[id_search_edge].Give_id_corner((dir1D)c);
                       if(id_point == c_SW)
                           quadrangles[id_quad].Set_id_interior_edge(SWdir2D,id_search_edge);
                       else if(id_point == c_SE)
                           quadrangles[id_quad].Set_id_interior_edge(SEdir2D,id_search_edge);
                       else if(id_point == c_NW)
                           quadrangles[id_quad].Set_id_interior_edge(NWdir2D,id_search_edge);
                       else if(id_point == c_NE)
                           quadrangles[id_quad].Set_id_interior_edge(NEdir2D,id_search_edge);
                     }
                  }
             }
          }
          else{
             // Exists exterior neighbor true
             quadrangles[id_quad].Set_exists_exterior(true);
             // Set exterior volume
             quadrangles[id_quad].Set_id_exterior_vol(num);
             quadrangles[id_quad].Set_typ_exterior_vol(hexahedronEl);

             // Set exterior faces
             // a) id of edges of Quad
             ed_W = quadrangles[id_quad].Give_id_edge(Wdir2D);
             ed_E = quadrangles[id_quad].Give_id_edge(Edir2D);
             ed_S = quadrangles[id_quad].Give_id_edge(Sdir2D);
             ed_N = quadrangles[id_quad].Give_id_edge(Ndir2D);
             // b) Iterate over all faces of hexahedron
             for(i=0;i<6;i++){
                 id_search_quad = hexahedra[num].Give_id_quadrangle((dir3D)i);
                 if(id_search_quad != id_quad){  // face is not the considered quad
                    for(ed = 0;ed<4;++ed){       // Iterate over all edges of faces
                                                  // and compare with edges of quad
                       id_edge = quadrangles[id_search_quad].Give_id_edge((dir2D)ed);
                       if(id_edge == ed_W)
                           quadrangles[id_quad].Set_id_exterior_face(Wdir2D,id_search_quad);
                       else if(id_edge == ed_E)
                           quadrangles[id_quad].Set_id_exterior_face(Edir2D,id_search_quad);
                       else if(id_edge == ed_S)
                           quadrangles[id_quad].Set_id_exterior_face(Sdir2D,id_search_quad);
                       else if(id_edge == ed_N)
                           quadrangles[id_quad].Set_id_exterior_face(Ndir2D,id_search_quad);
                     }
                  }
             }
             // c) Set type of exterior faces
             for(i = 0;i<4;++i){
                 quadrangles[id_quad].Set_typ_exterior_face((dir2D)i, quadrangleEl);
             }
             // Set exterior edges
             // a) id of corners of Quad
             c_SW = quadrangles[id_quad].Give_id_corner(SWdir2D);
             c_SE = quadrangles[id_quad].Give_id_corner(SEdir2D);
             c_NW = quadrangles[id_quad].Give_id_corner(NWdir2D);
             c_NE = quadrangles[id_quad].Give_id_corner(NEdir2D);
             // b) Iterate over all edges of hexahedron
             for(i=0;i<12;i++){
                 id_search_edge = hexahedra[num].Give_id_edge((Edges_cell)i);
                 if(id_search_edge != ed_W &&
                     id_search_edge != ed_E &&
                     id_search_edge != ed_S &&
                     id_search_edge != ed_N){  // edge belongs not to the considered quad
                     for(c = 0;c<2;++c){      // Iterate over all corners of edges
                                               // and compare with corners of quad
                       id_point = edges[id_search_edge].Give_id_corner((dir1D)c);
                       if(id_point == c_SW)
                           quadrangles[id_quad].Set_id_exterior_edge(SWdir2D,id_search_edge);
                       else if(id_point == c_SE)
                           quadrangles[id_quad].Set_id_exterior_edge(SEdir2D,id_search_edge);
                       else if(id_point == c_NW)
                           quadrangles[id_quad].Set_id_exterior_edge(NWdir2D,id_search_edge);
                       else if(id_point == c_NE)
                           quadrangles[id_quad].Set_id_exterior_edge(NEdir2D,id_search_edge);
                     }
                  }
             }
          }
      }
  }

  // construct neighbors of edges
  for(num=0;num<num_hexahedra;++num)
      for(i=0; i<12; ++i){
          id_edge = hexahedra[num].Give_id_edge((Edges_cell)i);
          edges[id_edge].Put_to_topf_neighbor((int)hexahedronEl, num);
      }
  for(num=0;num<num_quadrangles;++num)
      for(i=0; i<4; ++i){
          id_edge = quadrangles[num].Give_id_edge((dir2D)i);
          edges[id_edge].Put_to_topf_neighbor( (int)quadrangleEl, num);
      }
  for(id_edge=0;id_edge<num_edges;++id_edge){
      W_corner = edges[id_edge].Give_id_corner_W();
      E_corner = edges[id_edge].Give_id_corner_E();
      number_neighbors = edges[id_edge].Give_size_of_topf_neighbor();
      number_neighbors_hex = number_neighbors/2;
      edges[id_edge].Set_number_neighbors(number_neighbors);
      edges[id_edge].Set_number_neighbors_hex(number_neighbors_hex);
      edges[id_edge].Set_size_typ_M_el(number_neighbors);
      edges[id_edge].Set_size_typ_W_el(number_neighbors);
      edges[id_edge].Set_size_typ_E_el(number_neighbors);
      edges[id_edge].Set_size_id_M_el(number_neighbors);
      edges[id_edge].Set_size_id_W_el(number_neighbors);
      edges[id_edge].Set_size_id_E_el(number_neighbors);
      edges[id_edge].Set_size_loc_id_hex(number_neighbors_hex);
      edges[id_edge].Set_size_rel_neigh_locstiff(number_neighbors);
      for(num=0;num<number_neighbors;++num) {
          USPair = edges[id_edge].Find_in_topf_neighbor(num);
          edges[id_edge].Set_typ_M_el(num,(elementTyp)USPair.first);
          edges[id_edge].Set_id_M_el (num,USPair.second);
      }
      num_hex = 0;
      for(num=0;num<number_neighbors;++num) {
          neighbor_in_typ = edges[id_edge].Give_typ_M_el(num);
          neighbor_in_id = edges[id_edge].Give_id_M_el (num);
          if(neighbor_in_typ == hexahedronEl){
             edges[id_edge].Set_typ_W_el(num,quadrangleEl);
             edges[id_edge].Set_typ_E_el(num,quadrangleEl);
             edges[id_edge].Set_loc_id_hex(num_hex, NT, num);
             edges[id_edge].Set_rel_neigh_locstiff(num, num_hex);
             for(ed=0; ed<12; ++ed){
                  id_ed = hexahedra[neighbor_in_id].Give_id_edge((Edges_cell)ed);
                  if(id_ed == id_edge) ed_cell = (Edges_cell)ed;
             }
             if(hexahedra[neighbor_in_id].Give_id_corner(Transform(ed_cell, Ldir1D))<
             hexahedra[neighbor_in_id].Give_id_corner(Transform(ed_cell, Rdir1D))){
               edges[id_edge].Set_id_W_el
                (num,hexahedra[neighbor_in_id].Give_id_quadrangle(Trans_dir(ed_cell, Ldir1D)));
               edges[id_edge].Set_id_E_el
                (num,hexahedra[neighbor_in_id].Give_id_quadrangle(Trans_dir(ed_cell, Rdir1D)));
               id_global_ND = hexahedra[neighbor_in_id].Give_id_quadrangle(OrthoC(ed_cell));
               id_global_ST = hexahedra[neighbor_in_id].Give_id_quadrangle(OrthoD(ed_cell));
             }
             else{
               edges[id_edge].Set_id_W_el
                (num,hexahedra[neighbor_in_id].Give_id_quadrangle(Trans_dir(ed_cell, Rdir1D)));
               edges[id_edge].Set_id_E_el
                (num,hexahedra[neighbor_in_id].Give_id_quadrangle(Trans_dir(ed_cell, Ldir1D)));
               id_global_ND = hexahedra[neighbor_in_id].Give_id_quadrangle(OrthoD(ed_cell));
               id_global_ST = hexahedra[neighbor_in_id].Give_id_quadrangle(OrthoC(ed_cell));
             }
             for(num_q=0;num_q<number_neighbors;++num_q) {
                 if(edges[id_edge].Give_typ_M_el(num_q)==quadrangleEl){
                     if(edges[id_edge].Give_id_M_el(num_q)==id_global_ND)
                        edges[id_edge].Set_loc_id_hex(num_hex, ND, num_q);
                     else if(edges[id_edge].Give_id_M_el(num_q)==id_global_ST)
                        edges[id_edge].Set_loc_id_hex(num_hex, ST, num_q);
                 }
             }
             num_hex++;
             if(developer_version)
                if(num_hex>number_neighbors_hex)
                   cout << "error1 in construction_hexahedron_points_done!" << endl;
          }
          else if(neighbor_in_typ == quadrangleEl){
             edges[id_edge].Set_rel_neigh_locstiff(num, -1);
             edges[id_edge].Set_typ_W_el(num,edgeEl);
             edges[id_edge].Set_typ_E_el(num,edgeEl);
             for(e=0;e<4;++e){
                id_ed = quadrangles[neighbor_in_id].Give_id_edge((dir2D)e);
                if(id_ed != id_edge)
                   for(c=0;c<2;++c){
                       if(edges[id_ed].Give_id_corner((dir1D)c) == W_corner)
                          edges[id_edge].Set_id_W_el(num,id_ed);
                       else if(edges[id_ed].Give_id_corner((dir1D)c) == E_corner)
                          edges[id_edge].Set_id_E_el(num,id_ed);
                   }
             }
          }
      }
      if(developer_version)
         if(num_hex!=number_neighbors_hex)
            cout << "error2 in construction_hexahedron_points_done!" << endl;
  }

  //construcht neighbors of corners
  for(num=0;num<num_hexahedra;++num)
      for(i=0; i<8; ++i){
          id_point = hexahedra[num].Give_id_corner((dir3D_sons)i);
          points[id_point].Put_to_topf_neighbor((int)hexahedronEl, num);
      }
  for(num=0;num<num_quadrangles;++num)
      for(i=0; i<4; ++i){
          id_point = quadrangles[num].Give_id_corner((dir2D_sons)i);
          points[id_point].Put_to_topf_neighbor((int)quadrangleEl, num);
      }
  for(num=0;num<num_edges;++num)
      for(i=0; i<2; ++i){
          id_point = edges[num].Give_id_corner((dir1D)i);
          points[id_point].Put_to_topf_neighbor((int)edgeEl, num);
      }
  for(id_point=0;id_point<num_points;++id_point){
      number_neighbors = points[id_point].Give_size_of_topf_neighbor();
      number_neighbors_hex = 0;//vorlaeufig, Berechnung in for-schleife
      points[id_point].Set_number_neighbors(number_neighbors);
      points[id_point].Set_size_typ_el_neighbor(number_neighbors);
      points[id_point].Set_size_id_el_neighbor(number_neighbors);
      points[id_point].Set_size_rel_neigh_locstiff(number_neighbors);
      for(num=0;num<number_neighbors;++num) {
          USPair = points[id_point].Find_in_topf_neighbor(num);
          points[id_point].Set_typ_el_neighbor(num,(elementTyp)USPair.first);
          points[id_point].Set_id_el_neighbor (num,USPair.second);
          if((elementTyp)USPair.first == hexahedronEl) ++number_neighbors_hex;
      }
      points[id_point].Set_number_neighbors_hex(number_neighbors_hex);
      points[id_point].Set_size_loc_id_hex(number_neighbors_hex);
      num_hex = 0;
      for(num=0;num<number_neighbors;++num) {
          neighbor_in_typ = points[id_point].Give_typ_el_neighbor(num);
          neighbor_in_id = points[id_point].Give_id_el_neighbor(num);
          if(neighbor_in_typ == hexahedronEl){
             points[id_point].Set_loc_id_hex(num_hex, NET, num);
             points[id_point].Set_rel_neigh_locstiff(num, num_hex);
             for(c=0; c<8; ++c){
                  id_corner = hexahedra[neighbor_in_id].Give_id_corner((dir3D_sons)c);
                  if(id_corner == id_point) corner = (dir3D_sons)c;
             }
             Calculate_corresponding_corners(corner, SED_corner, NWD_corner, SWT_corner);
             Calculate_corresponding_edges(corner, SD_edge, WD_edge, SW_edge);
             id_global_SD = hexahedra[neighbor_in_id].Give_id_edge(SD_edge);
             id_global_WD = hexahedra[neighbor_in_id].Give_id_edge(WD_edge);
             id_global_SW = hexahedra[neighbor_in_id].Give_id_edge(SW_edge);
             id_global_W = hexahedra[neighbor_in_id].
                           Give_id_quadrangle(direction_from_to(SED_corner, corner));
             id_global_S = hexahedra[neighbor_in_id].
                           Give_id_quadrangle(direction_from_to(NWD_corner, corner));
             id_global_D = hexahedra[neighbor_in_id].
                           Give_id_quadrangle(direction_from_to(SWT_corner, corner));
             for(num_q=0;num_q<number_neighbors;++num_q) {
                 if(points[id_point].Give_typ_el_neighbor(num_q)==quadrangleEl){
                     if(points[id_point].Give_id_el_neighbor(num_q)==id_global_W)
                        points[id_point].Set_loc_id_hex(num_hex, NWT, num_q);
                     else if(points[id_point].Give_id_el_neighbor(num_q)==id_global_S)
                        points[id_point].Set_loc_id_hex(num_hex, SET, num_q);
                     else if(points[id_point].Give_id_el_neighbor(num_q)==id_global_D)
                        points[id_point].Set_loc_id_hex(num_hex, NED, num_q);
                 }
                 if(points[id_point].Give_typ_el_neighbor(num_q)==edgeEl){
                     if(points[id_point].Give_id_el_neighbor(num_q)==id_global_SD)
                        points[id_point].Set_loc_id_hex(num_hex, SED, num_q);
                     else if(points[id_point].Give_id_el_neighbor(num_q)==id_global_WD)
                        points[id_point].Set_loc_id_hex(num_hex, NWD, num_q);
                     else if(points[id_point].Give_id_el_neighbor(num_q)==id_global_SW)
                        points[id_point].Set_loc_id_hex(num_hex, SWT, num_q);
                 }
             }
             num_hex++;
             if(developer_version)
                if(num_hex>number_neighbors_hex)
                   cout << "error3 in construction_hexahedron_points_done!" << endl;
          }
          else points[id_point].Set_rel_neigh_locstiff(num, -1);
      }
      if(developer_version)
         if(num_hex!=number_neighbors_hex)
            cout << "error4 in construction_hexahedron_points_done!" << endl;
  }
  
  // set coordinates of quadrangles and edges
  Set_coordinates_quadrangles();
  Set_coordinates_edges();
}

// setze Koordinaten der Quadrangles von hexaeder ecken
void Unstructured_grid::Set_coordinates_quadrangles() {
  int      idpoint_hex[4];
  D3vector coord_hex[4];
  dir3D_sons co_hex;
  int id_quad;
  int idpoint_quad;
  
  for(int num=0;num<num_hexahedra;++num){
    for(int q=0; q<6; ++q) {
      for(int c=0;c<4;++c) {
        co_hex = corner_sur((dir3D)q,(dir2D_sons)c);
        idpoint_hex[c] = hexahedra[num].Give_id_corner((dir3D_sons)co_hex);
        coord_hex[c]   = hexahedra[num].Give_coord((dir3D_sons)co_hex);
      }
      id_quad = hexahedra[num].Give_id_quadrangle((dir3D)q);
      for(int cq=0;cq<4;++cq) {
        idpoint_quad = quadrangles[id_quad].Give_id_corner((dir2D_sons)cq);
        for(int c=0;c<4;++c) if(idpoint_hex[c]==idpoint_quad) 
          quadrangles[id_quad].Set_coord((dir2D_sons)cq,coord_hex[c]);
      }
    }
  }
}

// setze Koordinaten der Edges von hexaeder ecken
void Unstructured_grid::Set_coordinates_edges() {
  int      idpoint_hex[2];
  D3vector coord_hex[2];
  dir3D_sons co_hex;
  int id_edge;
  int idpoint_edge;
  
  for(int num=0;num<num_hexahedra;++num){
    for(int q=0; q<12; ++q) {
      for(int c=0;c<2;++c) {
        co_hex = Transform((Edges_cell)q,(dir1D)c);
        idpoint_hex[c] = hexahedra[num].Give_id_corner((dir3D_sons)co_hex);
        coord_hex[c]   = hexahedra[num].Give_coord((dir3D_sons)co_hex);
      }
      id_edge = hexahedra[num].Give_id_edge((Edges_cell)q);
      for(int cq=0;cq<2;++cq) {
        idpoint_edge = edges[id_edge].Give_id_corner((dir1D)cq);
        for(int c=0;c<2;++c) if(idpoint_hex[c]==idpoint_edge) 
          edges[id_edge].Set_coord((dir1D)cq,coord_hex[c]);
      }
    }
  }
}

void Unstructured_grid::construction_done() {
  int  color, num;
  int id_first_colored_edge;

  if(developer_version)
     if(construction_hexahedron_points_done_yn!=true)
        cout << " error in Unstructured_grid::construction_done !" << endl;

  // I. color edges with respect to degree of freedom
  color = 0;
  for(id_first_colored_edge=-2;id_first_colored_edge!=-1;++color) {
    id_first_colored_edge=-1;
    // 1. Find not colored edge
    for(num=0;num<num_edges && id_first_colored_edge==-1;++num) {
      if(edges[num].Give_color_of_edge()==-1) {
	edges[num].Set_color_of_edge(color);
	id_first_colored_edge = num;
      }
    }

    // 2. find edges with same color
    if(id_first_colored_edge!=-1)
      Set_recursive_color(id_first_colored_edge,color);
  }
  number_degree_of_freedom = color-1;

  // II. set transformation edges in hexahedra
  for(int id_hex=0;id_hex<num_hexahedra;++id_hex) {
    for(int i=0;i<12;++i) if(hexahedra[id_hex].transform[i] == NULL) {
        //todo : add if face transform at corresponding face is defined --> give_coord_edge() instead of transform_nothing ???

      // set transformation function
      hexahedra[id_hex].transform[i] =
	      edges[hexahedra[id_hex].Give_id_edge((Edges_cell)i)].transform;
      hexahedra[id_hex].setShiftPointer(edges[hexahedra[id_hex].Give_id_edge((Edges_cell)i)].getShiftPointer());
    };
  }

  // III. set all points
  all_points = new Marker(this);

  // IV. Zordering
  zordering = new Zordering(this);
  //zordering->Print(this);

  // V. Check, whether FACE or EDGE transform is used to find coordinates. Default is FALSE
  transformFromQuadrangle = true;

  //transformFromQuadrangle
  for (int qu = 0; qu < Give_number_quadrangles(); ++qu)
  {
      if( Give_quadrangle(qu)->transform == NULL)
      {
          transformFromQuadrangle = false;
      }
  }
  if (transformFromQuadrangle)
  { edge_to_quad_id.resize(Give_number_edges());
      edge_to_quad_dir.resize(Give_number_edges());
      edge_to_hex_id.resize(Give_number_edges());
  for (int iter = 0 ; iter <edge_to_quad_id.size(); iter++ )
  {
      edge_to_quad_id.at(iter) = -1;
     edge_to_hex_id.at(iter) = -1;
  }}

}

void Unstructured_grid::Set_recursive_color(int id_edge, int color) {
  int id_hex;
  dir3D main_dir;
  Edges_cell ec;
  int id_edge_next;

  // iterate neighbors of edge
  for(int i=0; i<edges[id_edge].Give_number_neighbors() ;++i) {
    // find volume element - missing prism
    if(edges[id_edge].Give_typ_M_el(i)==hexahedronEl) {
      id_hex = edges[id_edge].Give_id_M_el(i);
      // iterate edges of hexahedron with edge id_edge
      for(int ed=0;ed<12;++ed) {
	// find colored edge
	if(hexahedra[id_hex].Give_id_edge((Edges_cell)ed)==id_edge) {
	  main_dir = Trans_dir((Edges_cell)ed,Ldir1D);
          // iterate edges with the same color
	  for(int j=0;j<4;++j) {
	    ec = CalcEdge_cells(main_dir,(dir2D_sons)j);
	    if(ec!=(Edges_cell)ed) {
	      id_edge_next = hexahedra[id_hex].Give_id_edge(ec);
	      if(edges[id_edge_next].Give_color_of_edge()==-1) {
                // set color and recursive call
		edges[id_edge_next].Set_color_of_edge(color);
		Set_recursive_color(id_edge_next,color);
	      }
	    }
	  }
	}
      }
    }
  }
}

void Unstructured_grid::Set_transformation_edge(int idL, int idR,
						D3vector (*Transform)(double, double*),
						int shift_of_pointer) {
  // transform must be zero vector at 0,1
  // transform: [0,1] ->  D3vector
  if(developer_version){
     if(construction_hexahedron_points_done_yn!=true &&
        construction_done_yn==true)
        cout << " error 1 in Unstructured_grid::Set_transformation_edge !" << endl;
     if(idL>=idR)
        cout << " error 2 in Unstructured_grid::Set_transformation_edge !" << endl;
  }
  Edge_el* poiE = &edges[topf_edges.Give_id(idL,idR)];
  poiE->transform           = Transform;
  
  poiE->setShiftPointer(shift_of_pointer);
}

void Unstructured_grid::Set_transformation_edge_from_face(int idL, int idR, D3vector (*transform)(double, double, double *), int shift_of_pointer)
{

}


void Unstructured_grid::Set_transformation_edge_of_hexahedron(int id_hex,
							      Edges_cell ed,
						D3vector (*Transform)(double, double*)) {
  // transform must be zero vector at 0,1
  // transform: [0,1] ->  D3vector
  if(developer_version){
     if(construction_hexahedron_points_done_yn!=true &&
        construction_done_yn==true)
        cout << " error 1 in Unstructured_grid::Set_transformation_edge_of.. !" << endl;
     if(id_hex < 0 || id_hex >= num_hexahedra)
        cout << " error 2 in Unstructured_grid::Set_transformation_edge_of.. !" << endl;
  }
  hexahedra[id_hex].transform[(int)ed]   = Transform;

  assert(false); // das wurde noch nie angwendet. Es geht vielleicht hier aber nicht bei compose_grid.h
}

void Unstructured_grid::Set_transformation_face(int idSW, int idSE,
						int idNW, int idNE,
                        D3vector (*Transform)(double,double, double*),
                                                int shift_of_pointer) {
  // transform must be zero at (0,0)
  if(developer_version)
     if(construction_hexahedron_points_done_yn!=true &&
        construction_done_yn==true)
        cout << " error in Unstructured_grid::Set_transformation_face !" << endl;
 // cout << " Unstructured_grid::Set_transformation_face is being implemented ! ..." << endl;
  Quadrangle_el* quadEL = &quadrangles[topf_quadrangles.Give_id(idSW,idSE,idNW,idNE)];

  //cout << "id of quadrangle transform: " << topf_quadrangles.Give_id(idSW,idSE,idNW,idNE) << endl;

  quadEL->transform           = Transform;
  quadEL->setShiftPointer(shift_of_pointer);

  //Set_transformation_edge_from_face();

}

void Unstructured_grid::Set_transformation_hex(int id_hex,
                                                D3vector (*transform_h)(int ,double,double,double)) {
  if(developer_version)
    if(construction_hexahedron_points_done_yn!=true && construction_done_yn==true)
    cout << " error in Unstructured_grid::Set_transformation_hex !" << endl;
    hexahedra[id_hex].transform_hex = transform_h;

  assert(false); // das wurde noch nie angwendet. Es geht vielleicht hier aber nicht bei compose_grid.h  
}

void Unstructured_grid::Print() {
  cout << "Print ug 3D " << endl;
  cout << "---------------------------" << endl;
   for(int i=0; i<num_points;++i) {
     cout << "punkt: " << i << " : ";
     points[i].Give_coordinate().Print();
     cout << endl;
   }
}