meso2d.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:            #     #  ######   #####    #####    #####   ######
//                  ##   ##  #       #     #  #     #  #     #  #     #
//                  # # # #  #       #        #     #        #  #     #
//                  #  #  #  #####    #####   #     #   #####   #     #
//                  #     #  #             #  #     #  #        #     #
//                  #     #  #       #     #  #     #  #        #     #
//                  #     #  ######   #####    #####   #######  ######
//
// name:           meso2d.h
// Copyright:      Daniel Rypl
//                 Czech Technical University in Prague,
//                 Faculty of Civil Engineering, Department of Structural Mechanics, 
//                 Thakurova 7, 166 29 Prague, Czech Republic,
//                 email: drypl@fsv.cvut.cz
//
// language:       C, C++
// license:        This program is free software; you can redistribute it and/or modify
//                 it under the terms of the GNU Lesser General Public License as published by
//                 the Free Software Foundation; either version 2 of the License, or
//                 (at your option) any later version.
//
//                 This program is distributed in the hope that it will be useful,
//                 but WITHOUT ANY WARRANTY; without even the implied warranty of
//                 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//                 GNU Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************

#include <stdio.h>
#include <stdlib.h>
//#include <string.h>
//#include <ctype.h>
#include <math.h>
#include <stdarg.h>

#include "meso2d.h"

namespace meso2d {

//#ifndef ELIXIR
//#define ELIXIR
//#endif
//
//#ifdef ELIXIR
//#include "Esimple.h"
//#undef PI
//#undef NO
//#undef YES
//#endif


//    #define BUFFER_SIZE           1024
//    
//    #define REMARK               '#'
//    #define NEWLINE              '\n'
//    #define SPACE                ' '
//    #define TABULATOR            '\t'
//    

#define NO_ERROR             -1
#define WARNING               0
#define GENERAL_ERROR         1

//    #define FILE_OPEN_ERROR       2
//    #define FILE_READ_ERROR       3
//    #define FILE_WRITE_ERROR      4
//    #define MEMORY_ERROR          5
#define INPUT_ERROR           6
//    
//    

#define EPSILON               1.0e-10
#define EPSILON_T             1.0e-12

//    
//    #define ANGLE                 0.001         /* cos(89.94 degs) */
//    
//    #define PI                    M_PI          /* 3.14159265358979323846 */
//    
//    
//    
#define B_0(T) (1.0-2.0*(T)+(T)*(T))
#define B_1(T) (2.0*(1.0-(T))*(T))
#define B_2(T) ((T)*(T))
//    
#define dB_0(T) (-2.0+2.0*(T))
#define dB_1(T) (2.0-4.0*(T))
#define dB_2(T) (2.0*(T))
//    
//    
//    #define copy_vec(DEST, SRC)                {(DEST).x = (SRC).x; (DEST).y = (SRC).y;}
#define size_vec(VEC)                      ((VEC).x * (VEC).x + (VEC).y * (VEC).y)

#define dist_point(PNT1, PNT2)             (((PNT1).x - (PNT2).x) * ((PNT1).x - (PNT2).x) + \
                        ((PNT1).y - (PNT2).y) * ((PNT1).y - (PNT2).y))

//    #define aver_point(RES, PNT1, PNT2)        {(RES).x = ((PNT1).x + (PNT2).x) / 2.0; (RES).y = ((PNT1).y + (PNT2).y) / 2.0;}
//    
//    #define mul_vec(VEC, VALUE)                {(VEC).x *= (VALUE); (VEC).y *= (VALUE);}
//    #define div_vec(VEC, VALUE)                {(VEC).x /= (VALUE); (VEC).y /= (VALUE);}
//    #define add_vec(RES, VEC1, VEC2)           {(RES).x = (VEC1).x + (VEC2).x; (RES).y = (VEC1).y + (VEC2).y;}
#define sub_vec(RES, VEC1, VEC2)           {(RES).x = (VEC1).x - (VEC2).x; (RES).y = (VEC1).y - (VEC2).y;}
//    
//    
#define dot_product(VEC1, VEC2)            ((VEC1).x * (VEC2).x + (VEC1).y * (VEC2).y)
//    
//    
//    

//typedef enum logic logic;

//    
//    typedef struct point_rec point_rec;
//    
//    typedef struct ellipse_rec ellipse_rec;
//    typedef struct cylinder_rec cylinder_rec;
//    typedef struct crack_rec crack_rec;
//    
//    
//

enum logic {
  NO = 0,
  YES = 1
};

//    
//    

//    
//    

//    
//    
//    struct cylinder_rec{
//      point_rec center;
//      double radius, height;
//      double internal_msz, boundary_msz;
//      int property;
//      point_rec lower_left, lower_right, upper_left, upper_right;
//    };
//    
//    
//    struct crack_rec{
//      int id;
//      ellipse_rec *ellipse;
//      point_rec center, tip;
//      int interface, bc;
//      double length, angle, radius;
//      double crack_msz, tip_msz;
//      double cos_angle, sin_angle;
//      int quad;
//      int curve;
//      char bc_spec[128];
//    };
//    
//    
//    
//    static char *geom_file_name = NULL;
//    static char *t3d_in_file_name = NULL;
//    static char *ctrl_file_name = NULL;
//    
//    static FILE *geom_file = NULL;
//    static FILE *t3d_in_file = NULL;
//    static FILE *ctrl_file = NULL;
//    

static double epsilon = EPSILON;

//    
//    static double min_x, max_x;
//    static double min_y, max_y;
//    
//    static double gmsz = 0.0;
//    static int prop = 0;
//    
//    static char line_buffer[BUFFER_SIZE];
//    
//    static int vertex_id = 0;
//    static int curve_id = 0;
//    static int patch_id = 0;
//    
//    static int crack_id = 0;
//    
//    static logic default_msz = NO;
//    static logic bc = NO;
//    
//    
//    #ifdef ELIXIR
//    
//    #define FONT                  "-adobe-courier-bold-r-normal--*-120-*-*-m-*-iso8859-1"
//    
//    static EPixel green, red, yellow, blue, magenta, black;
//    static Font font;
//    
//    static void
//    draw_point(point_rec *pnt, EPixel color);
//    static void
//    draw_line(point_rec *pnt1, point_rec *pnt2, EPixel color);
//    static void
//    draw_elliptic_arc(point_rec *pnt1, point_rec *pnt2, point_rec *pnt3, double w1, double w2, double w3, EPixel color);
//    
//    static void
//    draw_cylinder(cylinder_rec *cylinder, EPixel color);
//    static void
//    draw_ellipse(ellipse_rec *ellipse, EPixel color);
//    static void
//    draw_crack(crack_rec *crack, EPixel color);
//    #endif
//    
//    
//    static void
//    write_ellipse(ellipse_rec *ellipse);
//    static void
//    write_crack(crack_rec *crack);
//    static void
//    write_cylinder(cylinder_rec *cylinder, ellipse_rec *ellipse_array, crack_rec *crack_array, int ellipses, int cracks);
//    
//    static logic
//    check_ellipse_inside_cylinder(ellipse_rec *ellipse, cylinder_rec *cylinder);

// Lukasovo - presunuto do headeru
// logic check_ellipse_ellipse_overlap (ellipse_rec *ellipse1, ellipse_rec *ellipse2);

//    static logic
//    check_crack_cylinder_intersection(crack_rec *crack, cylinder_rec *cylinder);
//    static logic
//    check_crack_ellipse_intersection(crack_rec *crack, ellipse_rec *ellipse);
//    static logic
//    check_crack_crack_intersection(crack_rec *crack1, crack_rec *crack2);
//    static logic
//    check_crack_inside_ellipse(crack_rec *crack, ellipse_rec *ellipse);
//    
//    static logic
//    check_ellipse_line_intersection(ellipse_rec *ellipse, point_rec *point1, point_rec *point2);

static logic check_point_inside_ellipse (point_rec *point, ellipse_rec *ellipse);

//    static logic
//    check_line_line_intersection(point_rec *pnta1, point_rec *pnta2, point_rec *pntb1, point_rec *pntb2);
//    

static int project_point_to_ellipse (point_rec *point, ellipse_rec *ellipse, point_rec *pnt, double *t_par);
static void get_ellipse_point_quadrant (ellipse_rec *ellipse, int quadrant, double t, point_rec *pnt);

//    static double
//    get_ellipse_gradient_quadrant(ellipse_rec *ellipse, int quadrant, double t, point_rec *grad);
//    static void
//    get_ellipse_normal_quadrant(ellipse_rec *ellipse, int quadrant, double t, point_rec *normal);
//    

static void get_ellipse_point (ellipse_rec *ellipse, double t, point_rec *pnt);
static double get_ellipse_gradient (ellipse_rec *ellipse, double t, point_rec *grad);

//    static void
//    get_ellipse_normal(ellipse_rec *ellipse, double t, point_rec *normal);
//    

static void transform_from_global_to_local (point_rec *global, point_rec *local, point_rec *origin, double cos_angle, double sin_angle);

static void transform_from_local_to_global (point_rec *global, point_rec *local, point_rec *origin, double cos_angle, double sin_angle);
//    
//    static char *
//    get_next_record(FILE *file, char *file_name, char *err_msg);
//    static char *
//    get_next_relevant_record(FILE *file, char *file_name, char *err_msg);
//    
//    static void
//    write_record(FILE *file, char *file_name, char *buffer);
//    

static void error_message (int exit_code, const char *format, ...);




void check_ellipse_rec_consistency (meso2d::ellipse_rec & L)
{
  if (L.max <= 0.0)          meso2d::error_message (INPUT_ERROR, "Positive major half-axis size of inclusion %d required", L.id + 1);
  if (L.min <= 0.0)          meso2d::error_message (INPUT_ERROR, "Positive minor half-axis size of inclusion %d required", L.id + 1);
  if (L.min > L.max)         meso2d::error_message (INPUT_ERROR, "Invalid half-axis ordering of inclusion %d", L.id + 1);

  if (L.internal_msz < 0.0)  meso2d::error_message (INPUT_ERROR, "Positive internal mesh size of inclusion %d required", L.id + 1);
  if (L.boundary_msz < 0.0)  meso2d::error_message (INPUT_ERROR, "Positive boundary mesh size of inclusion %d required", L.id + 1);
  if (L.property < 0)    meso2d::error_message (INPUT_ERROR, "Non-negative property of inclusion %d required", L.id + 1);
  if (L.interface < 0)   meso2d::error_message (INPUT_ERROR, "Non-negative interface of inclusion %d required", L.id + 1);
  if (L.bc < 0)          meso2d::error_message (INPUT_ERROR, "Non-negative boundary condition of inclusion %d required", L.id + 1);

  // Uhel je jiz v radianech
  if (L.angle < -6.28318531 || L.angle > 6.28318531)
    meso2d::error_message (INPUT_ERROR, "Invalid angle of inclusion %d", L.id + 1);

  if (L.property == 0 && L.interface != 0) {
    meso2d::error_message (WARNING, "Empty inclusion %d cannot have interface ==> interface ignored", L.id + 1);
    L.interface = 0;
  }
  if (L.interface == 0 && L.bc != 0) {
    meso2d::error_message (WARNING, "Inclusion without interface cannot be subjected to bc ==> bc ignored", L.id + 1);
    L.bc = 0;
  }

  return;
}

bool isZero (double zero, double a)
{
  if (a == 0.0)          return true;
  if (-zero<a && a<zero) return true;
  else                   return false;
}
void transfom_ellipse_rec (meso2d::ellipse_rec & L)
{
  meso2d::point_rec point;

  // Two ellipses with centroids on the x-axis both rotated 90/270 DEG with small distance between
  // are wrongly evaluated as overlaped.
  // I cannot find the reason, however the following rounding to 0.0 repairs this problem.
  L.cos_angle = cos (L.angle);
  L.sin_angle = sin (L.angle);

  double d = 0.0;
  double zer = 1.e-6;
  for (int i=0; i<3; i++) {
    if (i==0) d = -1.0;
    if (i==1) d =  0.0;
    if (i==2) d =  1.0;

    if (isZero(zer, L.cos_angle  - d))  L.cos_angle  = d;
    if (isZero(zer, L.sin_angle  - d))  L.sin_angle  = d;
  }


  L.patch = 0;

  point.x = L.max;
  point.y = 0.0;
  meso2d::transform_from_local_to_global (&(L.pnt0), &point, &(L.center), L.cos_angle, L.sin_angle);

  point.x = L.max;
  point.y = L.min;
  meso2d::transform_from_local_to_global (&(L.pnt1), &point, &(L.center), L.cos_angle, L.sin_angle);

  point.x = 0.0;
  point.y = L.min;
  meso2d::transform_from_local_to_global (&(L.pnt2), &point, &(L.center), L.cos_angle, L.sin_angle);

  return;
}


//    
//    int 
//    main(int argc, char** argv)
//    {
//      cylinder_rec cylinder;
//      ellipse_rec *ellipse_array = NULL, *ellipse = NULL, *el = NULL;
//      crack_rec *crack_array = NULL, *crack = NULL, *cr = NULL;
//      point_rec point, normal;
//      double t, product, msz;
//      int i, j, ellipse_id, ellipses, cracks, points, vertices;
//      char buffer[1024], text[1024];
//      
//    #ifdef ELIXIR
//      EView *view;
//      BOOLEAN success;
//    #endif
//    
//      if(argc < 4){
//        fprintf(stderr, "\n"); 
//        fprintf(stderr, "Usage: %s geom_file t3d_file ctrl_file [epsilon]\n\n", argv[0]);
//        fprintf(stderr, "       geom_file: description of geometry - input file\n");
//        fprintf(stderr, "       t3d_file : description of t3d model - output file\n");
//        fprintf(stderr, "       ctrl_file: t3d2merlin ctrl file (fragment) - output file\n\n");
//        fprintf(stderr, "       center radius height msz bmsz material\n");
//        fprintf(stderr, "       # of inclusions\n");
//        fprintf(stderr, "       center angle max min msz bmsz material interface bc\n");
//        fprintf(stderr, "       # of cracks\n");
//        fprintf(stderr, "       inclusion center angle length msz tipmsz interface bc radius\n"); 
//        fprintf(stderr, "       # of ctrl points\n");
//        fprintf(stderr, "       point msz\n\n");
//        exit(1);
//      }
//    
//      geom_file_name = argv[1];
//      t3d_in_file_name = argv[2];
//      ctrl_file_name = argv[3];
//    
//      if(argc > 4)epsilon = atof(argv[4]);
//      if(epsilon < 0.0){
//        error_message(WARNING, "Positive epsilon required ==> zero used by default");
//        epsilon = 0.0;
//      }
//      
//      if((geom_file = fopen(geom_file_name, "r")) == NULL)error_message(FILE_OPEN_ERROR, "File %s opening error", geom_file_name);
//      if((t3d_in_file = fopen(t3d_in_file_name, "w")) == NULL)error_message(FILE_OPEN_ERROR, "File %s opening error", t3d_in_file_name);
//      if((ctrl_file = fopen(ctrl_file_name, "w")) == NULL)error_message(FILE_OPEN_ERROR, "File %s opening error", ctrl_file_name);
//      
//    #ifdef ELIXIR
//      if(ESIBuildInterface(ESI_GRAPHIC_EDITOR_MASK, 1, argv) != 0){
//        fprintf(stderr, "Graphic interface error\n");
//        exit(1);
//      }
//      ESIPopup();
//      view = ElixirNewView("Meso2D", "Meso2D", "midnightblue", "white", 400, 400);
//      EMAttachView(ESIModel(), view);
//      EVSetViewOrientation(view, VIEW_ORIENT_TOP);
//      EVShowAxes(view, NO);
//      ESIHandleCmd("render ambient 0.1");
//      EVSetRenderMode(view, NORMAL_RENDERING);
//      
//      green = ColorGetPixelFromString("green", &success);
//      red = ColorGetPixelFromString("red", &success);
//      yellow = ColorGetPixelFromString("yellow", &success);
//      blue = ColorGetPixelFromString("blue", &success);
//      magenta = ColorGetPixelFromString("magenta", &success);
//      black = ColorGetPixelFromString("black", &success);
//      
//      font = FontGetFontFromString(FONT, &success);
//      if(success == NO)font = FontDefaultFont();
//    #endif
//    
//      get_next_relevant_record(geom_file, geom_file_name, "Missing problem description");
//      strcpy(text, line_buffer);
//    
//      get_next_relevant_record(geom_file, geom_file_name, "Missing cylinder description");
//      if(sscanf(line_buffer, "%lf %lf %lf %lf %lf %lf %d", &(cylinder.center.x), &(cylinder.center.y), 
//          &(cylinder.radius), &(cylinder.height), &(cylinder.internal_msz), &(cylinder.boundary_msz), 
//          &(cylinder.property)) != 7)
//        error_message(INPUT_ERROR, "Invalid record - cylinder description");
//      
//      if(cylinder.radius <= 0)
//        error_message(INPUT_ERROR, "Positive cylinder radius required");
//      if(cylinder.height <= 0)
//        error_message(INPUT_ERROR, "Positive cylinder height required");
//      if(cylinder.internal_msz < 0.0)
//        error_message(INPUT_ERROR, "Positive internal cylinder mesh size required");
//      if(cylinder.boundary_msz < 0.0)
//        error_message(INPUT_ERROR, "Positive boundary cylinder mesh size required");
//      if(cylinder.property <= 0)
//        error_message(INPUT_ERROR, "Positive cylinder property required");
//    
//      min_x = cylinder.center.x - cylinder.radius;
//      max_x = cylinder.center.x + cylinder.radius;
//      min_y = cylinder.center.y;
//      max_y = cylinder.center.y + cylinder.height;
//    
//      cylinder.lower_left.x = cylinder.upper_left.x = min_x;
//      cylinder.lower_right.x = cylinder.upper_right.x = max_x;
//      cylinder.lower_left.y = cylinder.lower_right.y = min_y;
//      cylinder.upper_left.y = cylinder.upper_right.y = max_y;
//    
//      gmsz = cylinder.internal_msz;
//      prop = cylinder.property;
//    
//    #ifdef ELIXIR
//      draw_cylinder(&cylinder, green);
//      EVFitAllIntoView(view);
//    #endif
//    
//      get_next_relevant_record(geom_file, geom_file_name, "Missing number of inclusions");
//      if(sscanf(line_buffer, "%d", &ellipses) != 1)
//        error_message(INPUT_ERROR, "Invalid record - number of inclusions");
//    
//      if(ellipses < 0)error_message(INPUT_ERROR, "Non-negative number of inclusions required");
//      
//      if(ellipses != 0){
//        if((ellipse_array = calloc(ellipses, sizeof(ellipse_rec))) == NULL)
//          error_message(MEMORY_ERROR, "Memory allocation error");
//    
//        for(i = 0; i < ellipses; i++){
//          ellipse = &(ellipse_array[i]);
//          ellipse -> id = i + 1;
//    
//          get_next_relevant_record(geom_file, geom_file_name, "Missing inclusion description");
//          if(sscanf(line_buffer, "%lf %lf %lf %lf %lf %lf %lf %d %d %d", &(ellipse -> center.x), &(ellipse -> center.y), 
//          &(ellipse -> angle), &(ellipse -> max), &(ellipse -> min),
//          &(ellipse -> internal_msz), &(ellipse -> boundary_msz), 
//          &(ellipse -> property), &(ellipse -> interface), &(ellipse -> bc)) != 10)
//      error_message(INPUT_ERROR, "Invalid record - inclusion %d description", i + 1);
//              
//          if(ellipse -> max <= 0.0)
//      error_message(INPUT_ERROR, "Positive major half-axis size of inclusion %d required", i + 1);
//          if(ellipse -> min <= 0.0)
//      error_message(INPUT_ERROR, "Positive minor half-axis size of inclusion %d required", i + 1);
//          if(ellipse -> min > ellipse -> max)
//      error_message(INPUT_ERROR, "Invalid half-axis ordering of inclusion %d", i + 1);
//          if(ellipse -> internal_msz < 0.0)
//      error_message(INPUT_ERROR, "Positive internal mesh size of inclusion %d required", i + 1);
//          if(ellipse -> boundary_msz < 0.0)
//      error_message(INPUT_ERROR, "Positive boundary mesh size of inclusion %d required", i + 1);
//          if(ellipse -> property < 0)
//      error_message(INPUT_ERROR, "Non-negative property of inclusion %d required", i + 1);
//          if(ellipse -> interface < 0)
//      error_message(INPUT_ERROR, "Non-negative interface of inclusion %d required", i + 1);
//          if(ellipse -> bc < 0)
//      error_message(INPUT_ERROR, "Non-negative boundary condition of inclusion %d required", i + 1);
//          if(ellipse -> angle < -360.0 || ellipse -> angle > 360.0)
//      error_message(INPUT_ERROR, "Invalid angle of inclusion %d", i + 1);
//    
//          if(ellipse -> property == 0 && ellipse -> interface != 0){
//      error_message(WARNING, "Empty inclusion %d cannot have interface ==> interface ignored", i + 1);
//      ellipse -> interface = 0;
//          }
//          if(ellipse -> interface == 0 && ellipse -> bc != 0){
//      error_message(WARNING, "Inclusion without interface cannot be subjected to bc ==> bc ignored", i + 1);
//      ellipse -> bc = 0;
//          }
//              
//          ellipse -> angle = ellipse -> angle / 180.0 * PI;
//          ellipse -> cos_angle = cos(ellipse -> angle);
//          ellipse -> sin_angle = sin(ellipse -> angle);
//    
//          ellipse -> crack = NULL;
//          ellipse -> patch = 0;
//    
//          point.x = ellipse -> max;
//          point.y = 0.0;
//          transform_from_local_to_global(&(ellipse -> pnt0), &point, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//              
//          point.x = ellipse -> max;
//          point.y = ellipse -> min;
//          transform_from_local_to_global(&(ellipse -> pnt1), &point, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//              
//          point.x = 0.0;
//          point.y = ellipse -> min;
//          transform_from_local_to_global(&(ellipse -> pnt2), &point, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//    
//          if(check_ellipse_inside_cylinder(ellipse, &cylinder) == NO){
//      error_message(WARNING, "Inclusion %d intersects cylinder ==> inclusion %d ignored", i + 1, i + 1);
//      ellipse -> id = -ellipse -> id;
//    
//    #ifdef ELIXIR
//      draw_ellipse(ellipse, blue);
//    #endif
//    
//      continue;
//          }
//              
//          for(j = 0; j < i; j++){
//      el = &(ellipse_array[j]);
//      if(el -> id < 0)continue;
//      if(check_ellipse_ellipse_overlap(ellipse, el) == YES){
//        error_message(WARNING, "Inclusion %d intersects inclusion %d ==> inclusion %d ignored", i + 1, j + 1, i + 1);
//        ellipse -> id = -ellipse -> id;
//                      
//    #ifdef ELIXIR
//        draw_ellipse(ellipse, blue);
//    #endif
//    
//        break;
//      }
//          }
//          if(ellipse -> id < 0)continue;
//    
//    #ifdef ELIXIR
//          if(ellipse -> interface == 0)
//      draw_ellipse(ellipse, red);
//          else
//      draw_ellipse(ellipse, magenta);
//    #endif
//    
//        }
//      }
//    
//      get_next_relevant_record(geom_file, geom_file_name, "Missing number of cracks");
//      if(sscanf(line_buffer, "%d", &cracks) != 1)
//        error_message(INPUT_ERROR, "Invalid record - number of cracks");
//    
//      if(cracks < 0)error_message(INPUT_ERROR, "Non-negative number of cracks required");
//    
//      if(cracks != 0){
//        if((crack_array = calloc(cracks, sizeof(crack_rec))) == NULL)
//          error_message(MEMORY_ERROR, "Memory allocation error");
//    
//        for(i = 0; i < cracks; i++){
//          crack = &(crack_array[i]);
//          crack -> id = i + 1;
//    
//          get_next_relevant_record(geom_file, geom_file_name, "Missing crack description");
//          if(sscanf(line_buffer, "%d %lf %lf %lf %lf %lf %lf %d %d %lf", &ellipse_id, &(crack -> center.x), &(crack -> center.y), 
//          &(crack -> angle), &(crack -> length), &(crack -> crack_msz), &(crack -> tip_msz), 
//          &(crack -> interface), &(crack -> bc), &(crack -> radius)) != 10)
//      error_message(INPUT_ERROR, "Invalid record - crack %d description", i + 1);
//              
//          if(ellipse_id < 0 || ellipse_id > ellipses)
//      error_message(INPUT_ERROR, "Invalid inclusion number for crack %d", i + 1);
//          if(crack -> length <= 0.0)
//      error_message(INPUT_ERROR, "Positive length of crack %d required", i + 1);
//          if(crack -> crack_msz < 0.0)
//      error_message(INPUT_ERROR, "Positive mesh size of crack %d required", i + 1);
//          if(crack -> tip_msz < 0.0)
//      error_message(INPUT_ERROR, "Positive tip mesh size of crack %d required", i + 1);
//          if(crack -> angle < -360.0 || crack -> angle > 360.0)
//      error_message(INPUT_ERROR, "Invalid angle of crack %d", i + 1);
//          if(crack -> interface < 0)
//      error_message(INPUT_ERROR, "Non-negative interface of crack %d required", i + 1);
//          if(crack -> bc < 0)
//      error_message(INPUT_ERROR, "Non-negative boundary condition of crack %d required", i + 1);
//          if(crack -> radius < 0.0)
//      error_message(INPUT_ERROR, "Non-negative radius of crack %d required", i + 1);
//          /*          
//                  if(crack -> interface == 0 && crack -> radius == 0.0)
//                  error_message(WARNING, "Zero contour path radius assigned to crack %d without interface", i + 1);
//                  if(crack -> interface != 0 && crack -> radius != 0.0)
//                  error_message(WARNING, "Non-zero contour path radius assigned to crack %d with interface", i + 1);
//          */
//          crack -> angle = crack -> angle / 180.0 * PI;
//          crack -> cos_angle = cos(crack -> angle);
//          crack -> sin_angle = sin(crack -> angle);
//    
//          crack -> ellipse = NULL;
//          crack -> curve = 0;
//    
//          if(ellipse_id == 0){
//      ellipse = NULL;
//    
//      crack -> length *= 2.0;
//    
//      if(crack -> radius > crack -> length * 0.25){
//        error_message(WARNING, "Contour radius of crack %d too large ==> radius ignored", i + 1);
//        crack -> radius = 0.0;
//      }
//    
//      crack -> center.x = crack -> center.x - crack -> cos_angle * crack -> length * 0.5;
//      crack -> center.y = crack -> center.y - crack -> sin_angle * crack -> length * 0.5;
//    
//      crack -> tip.x = crack -> center.x + crack -> cos_angle * crack -> length;
//      crack -> tip.y = crack -> center.y + crack -> sin_angle * crack -> length;
//    
//      for(j = 0; j < ellipses; j++){
//        el = &(ellipse_array[j]);
//        if(el -> id < 0)continue;
//        if(check_crack_ellipse_intersection(crack, el) == YES){
//          error_message(WARNING, "Crack %d intersects inclusion %d ==> crack %d ignored", i + 1, j + 1, i + 1);
//          crack -> id = -crack -> id;
//    
//    #ifdef ELIXIR
//          draw_crack(crack, blue);
//    #endif
//    
//          break;
//        }
//        if(check_crack_inside_ellipse(crack, el) == YES){
//          error_message(WARNING, "Crack %d inside inclusion %d ==> crack %d ignored", i + 1, j + 1, i + 1);
//          crack -> id = -crack -> id;
//    
//    #ifdef ELIXIR
//          draw_crack(crack, blue);
//    #endif
//    
//          break;
//        }
//      }
//      if(crack -> id < 0)continue;
//          }
//          else{
//      ellipse = &(ellipse_array[ellipse_id - 1]);
//    
//      if(crack -> radius > crack -> length * 0.75){
//        error_message(WARNING, "Contour radius of crack %d too large ==> radius ignored", i + 1);
//        crack -> radius = 0.0;
//      }
//    
//      crack -> quad = project_point_to_ellipse(&(crack -> center), ellipse, &point, &t);
//      copy_vec(crack -> center, point);
//      get_ellipse_normal_quadrant(ellipse, crack -> quad, t, &normal);
//    
//      product = crack -> cos_angle * normal.x + crack -> sin_angle * normal.y;
//      if(fabs(product) < ANGLE){
//        error_message(WARNING, "Crack %d is tangent to inclusion %d ==> crack %d ignored", i + 1, ellipse_id, i + 1);
//        crack -> id = -crack -> id;
//    
//    #ifdef ELIXIR
//        crack -> tip.x = crack -> center.x + crack -> cos_angle * crack -> length;
//        crack -> tip.y = crack -> center.y + crack -> sin_angle * crack -> length;
//        draw_crack(crack, blue);
//    #endif
//    
//        continue;
//      }
//    
//      if(product < 0.0){
//        /*
//          error_message(WARNING, "Crack %d intersects inclusion %d ==> crack %d ignored", i + 1, ellipse_id, i + 1);
//          crack -> id = -crack -> id;
//    
//          #ifdef ELIXIR
//          crack -> tip.x = crack -> center.x + crack -> cos_angle * crack -> length;
//          crack -> tip.y = crack -> center.y + crack -> sin_angle * crack -> length;
//          draw_crack(crack, blue);
//          #endif
//    
//          continue;
//        */
//        error_message(WARNING, "Crack %d points inside inclusion %d ==> crack %d reverted", i + 1, ellipse_id, i + 1);
//    
//        crack -> angle += PI;
//        if(crack -> angle > 2.0 * PI)crack -> angle -= 2.0 * PI;
//        crack -> cos_angle *= -1.0;
//        crack -> sin_angle *= -1.0;
//      }
//                  
//      crack -> tip.x = crack -> center.x + crack -> cos_angle * crack -> length;
//      crack -> tip.y = crack -> center.y + crack -> sin_angle * crack -> length;
//    
//      if(ellipse -> id < 0){
//        error_message(WARNING, "Crack %d is on rejected inclusion ==> crack %d ignored", i + 1, i + 1);
//        crack -> id = -crack -> id;
//    
//    #ifdef ELIXIR
//        draw_crack(crack, blue);
//    #endif
//    
//        continue;
//      }
//    
//      if(ellipse -> crack != NULL){
//        if(ellipse -> crack -> id > 0){
//          error_message(WARNING, "Inclusion %d associated with crack %d ==> crack %d ignored", 
//                ellipse -> id, ellipse -> crack -> id, i + 1);
//          crack -> id = -crack -> id;
//    
//    #ifdef ELIXIR
//          draw_crack(crack, blue);
//    #endif
//    
//          continue;
//        }
//      }
//          }
//    
//          if(check_crack_cylinder_intersection(crack, &cylinder) == YES){
//      error_message(WARNING, "Crack %d intersects cylinder ==> crack %d ignored", i + 1, i + 1);
//      crack -> id = -crack -> id;
//    
//    #ifdef ELIXIR
//      draw_crack(crack, blue);
//    #endif
//    
//      continue;
//          }
//    
//          for(j = 0; j < i; j++){
//      cr = &(crack_array[j]);
//      if(cr -> id < 0)continue;
//      if(check_crack_crack_intersection(crack, cr) == YES){
//        error_message(WARNING, "Crack %d intersects crack %d ==> crack %d ignored", i + 1, j + 1, i + 1);
//        crack -> id = -crack -> id;
//      
//    #ifdef ELIXIR
//        draw_crack(crack, blue);
//    #endif
//    
//        break;
//      }
//          }
//          if(crack -> id < 0)continue;
//    
//          for(j = 0; j < ellipses; j++){
//      el = &(ellipse_array[j]);
//      if(el -> id < 0)continue;
//      if(el == ellipse)continue;
//      if(check_crack_ellipse_intersection(crack, el) == YES){
//        error_message(WARNING, "Crack %d intersects inclusion %d ==> crack %d ignored", i + 1, j + 1, i + 1);
//        crack -> id = -crack -> id;
//    
//    #ifdef ELIXIR
//        draw_crack(crack, blue);
//    #endif
//    
//        break;
//      }
//          }
//          if(crack -> id < 0)continue;
//              
//    #ifdef ELIXIR
//          draw_crack(crack, yellow);
//    #endif
//    
//          if(ellipse != NULL){
//      crack -> ellipse = ellipse;
//      ellipse -> crack = crack;
//          }
//        }
//      }
//    
//      /* note: text contains \n */
//      sprintf(buffer, "#\n# %s#\n\n", text);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      sprintf(buffer, "StartMesh\n");
//      write_record(ctrl_file, ctrl_file_name, buffer);
//    
//      /* note: text contains \n */
//      sprintf(buffer, "\n#\n# %s#\n\n", text);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//    
//      vertex_id = 10;
//      curve_id = 10;
//      patch_id = 10;
//    
//      for(j = 0; j < ellipses; j++){
//        ellipse = &(ellipse_array[j]);
//        write_ellipse(ellipse);
//      }
//    
//      for(j = 0; j < cracks; j++){
//        crack = &(crack_array[j]);
//        write_crack(crack);
//      }
//    
//      vertices = vertex_id;
//    
//      vertex_id = 0;
//      curve_id = 0;
//      patch_id = 0;
//    
//      write_cylinder(&cylinder, ellipse_array, crack_array, ellipses, cracks);
//    
//      sprintf(buffer, "EndMesh\n");
//      write_record(ctrl_file, ctrl_file_name, buffer);
//    
//      if(bc == YES){
//        sprintf(buffer, "\n#StartIncrement\n");
//        write_record(ctrl_file, ctrl_file_name, buffer);
//    
//        for(j = 0; j < ellipses; j++){
//          ellipse = &(ellipse_array[j]);
//          if(ellipse -> id < 0)continue;
//          if(ellipse -> bc != 0){
//      sprintf(buffer, "\n## inclusion no %d (inclusion and matrix surface)\n\n", ellipse -> id);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//      write_record(ctrl_file, ctrl_file_name, ellipse -> bc_ellipse_spec);
//      write_record(ctrl_file, ctrl_file_name, ellipse -> bc_matrix_spec);
//      sprintf(buffer, "#Tractions 0.0 0.0\n");
//      write_record(ctrl_file, ctrl_file_name, buffer);
//          }
//        }
//    
//        for(j = 0; j < cracks; j++){
//          crack = &(crack_array[j]);
//          if(crack -> id < 0)continue;
//          if(crack -> bc != 0){
//      sprintf(buffer, "\n## crack no %d\n\n", crack -> id);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//      write_record(ctrl_file, ctrl_file_name, crack -> bc_spec);
//      sprintf(buffer, "#Tractions 0.0 0.0\n");
//      write_record(ctrl_file, ctrl_file_name, buffer);
//          }
//        }
//    
//        sprintf(buffer, "\n#EndIncrement\n");
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//    
//      get_next_relevant_record(geom_file, geom_file_name, "Missing number of ctrl  points");
//      if(sscanf(line_buffer, "%d", &points) != 1)
//        error_message(INPUT_ERROR, "Invalid record - number of ctrl points");
//    
//      if(points < 0)error_message(INPUT_ERROR, "Non-negative number of ctrl points required");
//      
//      if(points != 0){
//        vertex_id = vertices;
//    
//        sprintf(buffer, "\n\n# control points\n\n");
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//        for(i = 0; i < points; i++){
//          get_next_relevant_record(geom_file, geom_file_name, "Missing ctrl point description");
//          if(sscanf(line_buffer, "%lf %lf %lf", &(point.x), &(point.y), &msz) != 3)
//      error_message(INPUT_ERROR, "Invalid record - ctrl point %d description", i + 1);
//    
//          if(msz < 0.0)
//      error_message(INPUT_ERROR, "Positive mesh size of ctrl point %d required", i + 1);
//              
//    #ifdef ELIXIR
//          draw_point(&point, green);
//    #endif
//    
//          vertex_id++;
//          if(msz == 0.0)
//      sprintf(buffer, "vertex %d virtual xyz %e %e %e\n", vertex_id, point.x, point.y, 0.0);
//          else
//      sprintf(buffer, "vertex %d virtual xyz %e %e %e size %e\n", vertex_id, point.x, point.y, 0.0, msz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//        }
//      }
//    
//      if(ellipse_array != NULL)free(ellipse_array);
//      if(crack_array != NULL)free(crack_array);
//    
//      if(default_msz == YES){
//        /*
//          error_message(WARNING, "Default mesh size (option -d) must be used");
//          error_message(NO_ERROR, "         when running t3d with %s", t3d_in_file_name);
//        */
//        sprintf(buffer, "\n\n#\n# t3d -i %s -o %s.out -p 24 -d ???\n#\n\n", t3d_in_file_name, t3d_in_file_name);
//      }
//      else
//        sprintf(buffer, "\n\n#\n# t3d -i %s -o %s.out -p 24\n#\n\n", t3d_in_file_name, t3d_in_file_name);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      fclose(geom_file);
//      fclose(t3d_in_file);
//      fclose(ctrl_file);
//    
//    #ifdef ELIXIR
//      ESIEventLoop(YES, "Completed");
//    #endif
//    
//      return(0);
//    }
//    
//    
//    
//    static void
//    write_ellipse(ellipse_rec *ellipse)
//    {
//      crack_rec *crack = NULL;
//      point_rec glob_point[4], loc_point[4], glob_pnt[4], loc_pnt[4], glob_p[2], loc_p[2];
//      point_rec *g_point[4], *g_pnt[4];
//      double weight_q = 0.707106781, weight_h = 0.3333333333, bmsz, imsz, cmsz, tmsz, min_msz, dist1, dist2, dist;
//      int i, k, quadrant, ver_id1, ver_id2, vertices, curves, crack_cur_id;
//      int matrix_crack_ver_id, interface_crack_ver_id, virtual_curve_id, tip_crack_ver_id, orig_matrix_crack_ver_id;
//      char buffer[1024], tmp_buffer[1024];
//    
//      if(ellipse -> id < 0){
//        sprintf(buffer, "# inclusion no %d (ignored)\n\n", -ellipse -> id);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//        return;
//      }
//    
//      if(ellipse -> property == 0)ellipse -> interface = 0;
//    
//      bmsz = ellipse -> boundary_msz;
//      imsz = ellipse -> internal_msz;
//    
//      /* vertices */
//      loc_point[0].x = ellipse -> max;
//      loc_point[0].y = 0.0;
//    
//      loc_point[1].x = 0.0;
//      loc_point[1].y = ellipse -> min;
//    
//      loc_point[2].x = -ellipse -> max;
//      loc_point[2].y = 0.0;
//    
//      loc_point[3].x = 0.0;
//      loc_point[3].y = -ellipse -> min;
//    
//      /* ctrl polygon */
//      loc_pnt[0].x = ellipse -> max;
//      loc_pnt[0].y = ellipse -> min;
//          
//      loc_pnt[1].x = -ellipse -> max;
//      loc_pnt[1].y = ellipse -> min;
//          
//      loc_pnt[2].x = -ellipse -> max;
//      loc_pnt[2].y = -ellipse -> min;
//          
//      loc_pnt[3].x = ellipse -> max;
//      loc_pnt[3].y = -ellipse -> min;
//          
//      for(i = 0; i < 4; i++){
//        transform_from_local_to_global(&(glob_point[i]), &(loc_point[i]), &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//        transform_from_local_to_global(&(glob_pnt[i]), &(loc_pnt[i]), &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//      }
//    
//      sprintf(buffer, "\n# inclusion no %d ", ellipse -> id);
//      if(ellipse -> property == 0)
//        strcat(buffer, "(hole)");
//      else{
//        if(ellipse -> interface == 0)
//          sprintf(tmp_buffer, "(material %d)", ellipse -> property);
//        else
//          sprintf(tmp_buffer, "(material %d, interface %d)", ellipse -> property, ellipse -> interface);
//        strcat(buffer, tmp_buffer);
//      }
//    
//      crack = ellipse -> crack;
//      if(crack != NULL){
//        if(crack -> interface == 0)
//          sprintf(tmp_buffer, " + crack no %d\n\n", crack -> id);
//        else
//          sprintf(tmp_buffer, " + crack no %d (interface %d)\n\n", crack -> id, crack -> interface);
//        strcat(buffer, tmp_buffer);
//      
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//    
//        cmsz = crack -> crack_msz;
//        tmsz = crack -> tip_msz;
//    
//        quadrant = crack -> quad;
//    
//        /* check whether crack is close to endpoints of quadrant */
//        dist1 = dist_point(crack -> center, glob_point[quadrant - 1]);
//        dist2 = dist_point(crack -> center, glob_point[quadrant % 4]);
//        dist = (dist1 < dist2) ? dist1 : dist2;
//    
//        min_msz = ellipse -> min / 2.0;
//        if(gmsz < min_msz && gmsz != 0.0)min_msz = gmsz;
//        if(imsz < min_msz && imsz != 0.0)min_msz = imsz;
//        if(bmsz < min_msz && bmsz != 0.0)min_msz = bmsz;
//        if(cmsz < min_msz && cmsz != 0.0)min_msz = cmsz;
//        min_msz *= 2.0;
//    
//        if(dist < min_msz * min_msz){
//    
//          /* setup the points that the half-ellipse is always in the first and second quadrants */
//          if(dist1 < dist2){
//      for(i = 0; i < 4; i++){
//        g_point[i] = &(glob_point[(quadrant - 2 + i + 4) % 4]);
//        g_pnt[i] = &(glob_pnt[(quadrant - 2 + i + 4) % 4]);
//      }
//    
//      /* calculate ctrl polygon of half-ellipse */
//      if(quadrant == 2){
//        loc_p[0].x = ellipse -> max;
//        loc_p[0].y = 2.0 * ellipse -> min;
//                      
//        loc_p[1].x = -ellipse -> max;
//        loc_p[1].y = 2.0 * ellipse -> min;
//      }
//      if(quadrant == 3){
//        loc_p[0].x = -2.0 * ellipse -> max;
//        loc_p[0].y = ellipse -> min;
//                      
//        loc_p[1].x = -2.0 * ellipse -> max;
//        loc_p[1].y = -ellipse -> min;
//      }
//      if(quadrant == 4){
//        loc_p[0].x = -ellipse -> max;
//        loc_p[0].y = -2.0 * ellipse -> min;
//                      
//        loc_p[1].x = ellipse -> max;
//        loc_p[1].y = -2.0 * ellipse -> min;
//      }
//      if(quadrant == 1){
//        loc_p[0].x = 2.0 * ellipse -> max;
//        loc_p[0].y = -ellipse -> min;
//                      
//        loc_p[1].x = 2.0 * ellipse -> max;
//        loc_p[1].y = ellipse -> min;
//      }
//          }
//          else{
//      for(i = 0; i < 4; i++){
//        g_point[i] = &(glob_point[(quadrant - 1 + i + 4) % 4]);
//        g_pnt[i] = &(glob_pnt[(quadrant - 1 + i + 4) % 4]);
//      }
//                  
//      /* calculate ctrl polygon of half-ellipse */
//      if(quadrant == 1){
//        loc_p[0].x = ellipse -> max;
//        loc_p[0].y = 2.0 * ellipse -> min;
//                      
//        loc_p[1].x = -ellipse -> max;
//        loc_p[1].y = 2.0 * ellipse -> min;
//      }
//      if(quadrant == 2){
//        loc_p[0].x = -2.0 * ellipse -> max;
//        loc_p[0].y = ellipse -> min;
//                      
//        loc_p[1].x = -2.0 * ellipse -> max;
//        loc_p[1].y = -ellipse -> min;
//      }
//      if(quadrant == 3){
//        loc_p[0].x = -ellipse -> max;
//        loc_p[0].y = -2.0 * ellipse -> min;
//                      
//        loc_p[1].x = ellipse -> max;
//        loc_p[1].y = -2.0 * ellipse -> min;
//      }
//      if(quadrant == 4){
//        loc_p[0].x = 2.0 * ellipse -> max;
//        loc_p[0].y = -ellipse -> min;
//                      
//        loc_p[1].x = 2.0 * ellipse -> max;
//        loc_p[1].y = ellipse -> min;
//      }
//          }
//    
//          transform_from_local_to_global(&(glob_p[0]), &(loc_p[0]), &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//          transform_from_local_to_global(&(glob_p[1]), &(loc_p[1]), &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//                  
//          vertices = 5;
//          curves = 4 + 1;   /* one is virtual */
//    
//          for(k = 0; k < 4; k++){
//    
//      /* k = 0 inclusion boundary
//         k = 1 matrix boundary
//         k = 2 inclusion interface
//         k = 3 matrix interface */
//    
//      if(k == 1){
//        sprintf(buffer, "MasterSlave");
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//      if(k == 2){
//        sprintf(buffer, "Interface %d", ellipse -> interface);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//    
//      for(i = 0; i < 4; i++){
//        if(i == 1)continue;
//        vertex_id++;
//        if(bmsz == 0.0){
//          default_msz = YES;
//          sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, g_point[i] -> x, g_point[i] -> y, 0.0);
//        }
//        else
//          sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, g_point[i] -> x, g_point[i] -> y, 0.0, bmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//        if(k == 1){
//          if(ellipse -> interface == 0)
//            sprintf(buffer, " Vertex %d %d", vertex_id - vertices, vertex_id);
//          else
//            sprintf(buffer, " Vertex %d %d Vertex %d %d", 
//                vertex_id - vertices, vertex_id + vertices, vertex_id, vertex_id + 2 * vertices);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//      }
//    
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      /* make virtual curve over the first two quadrants */
//      curve_id++;
//      sprintf(buffer, "curve %d virtual order 4 vertex %d %d\n", curve_id, vertex_id - 2, vertex_id - 1);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      virtual_curve_id = curve_id;
//                  
//      if(bmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "polygon 1 xyz %e %e %e weight %e\n", glob_p[0].x, glob_p[0].y, 0.0, weight_h);
//      }
//      else
//        sprintf(buffer, "polygon 1 xyz %e %e %e weight %e size %e\n", glob_p[0].x, glob_p[0].y, 0.0, weight_h, bmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      if(bmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "polygon 2 xyz %e %e %e weight %e\n", glob_p[1].x, glob_p[1].y, 0.0, weight_h);
//      }
//      else
//        sprintf(buffer, "polygon 2 xyz %e %e %e weight %e size %e\n", glob_p[1].x, glob_p[1].y, 0.0, weight_h, bmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      for(i = 2; i < 4; i++){
//        ver_id1 = vertex_id - 3 + i;
//        ver_id2 = ver_id1 + 1;
//        if(i == 3)ver_id2 -= 3;
//                  
//        curve_id++;
//        sprintf(buffer, "curve %d order 3 vertex %d %d\n", curve_id, ver_id1, ver_id2);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//                      
//        if(bmsz == 0.0){
//          default_msz = YES;
//          sprintf(buffer, "polygon 1 xyz %e %e %e weight %e\n", g_pnt[i] -> x, g_pnt[i] -> y, 0.0, weight_q);
//        }
//        else
//          sprintf(buffer, "polygon 1 xyz %e %e %e weight %e size %e\n", g_pnt[i] -> x, g_pnt[i] -> y, 0.0, weight_q, bmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//        if(k == 1){
//          if(ellipse -> interface == 0)
//            sprintf(buffer, " Curve %d %d", curve_id - curves, curve_id);
//          else
//            sprintf(buffer, " Curve %d %d Curve %d %d", 
//                curve_id - curves, curve_id + curves, curve_id, curve_id + 2 * curves);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//        if(k == 2){
//          sprintf(buffer, " Curve %d %d", curve_id + curves, curve_id);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//      }
//                  
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      /* split virtual curve */           
//      vertex_id++;
//      sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, 
//          virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      if(k == 0 || k == 1){
//        matrix_crack_ver_id = vertex_id;
//        interface_crack_ver_id = 0;
//      }
//      if(k == 3)interface_crack_ver_id = vertex_id;
//    
//      if(k == 1){
//        if(ellipse -> interface == 0)
//          sprintf(buffer, " Vertex %d %d", vertex_id - vertices, vertex_id);
//        else
//          sprintf(buffer, " Vertex %d %d Vertex %d %d", 
//              vertex_id - vertices, vertex_id + vertices, vertex_id, vertex_id + 2 * vertices);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//    
//      vertex_id++;
//      if(k == 1 || k == 3 || ellipse -> property == 0){
//        sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, 
//            virtual_curve_id);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//        if(k == 1){
//          if(ellipse -> interface == 0)
//            sprintf(buffer, " Vertex %d %d", vertex_id - 1 - vertices, vertex_id);
//          else{
//            if(ellipse -> property == 0)
//          sprintf(buffer, " Vertex %d %d Vertex %d %d",
//              vertex_id - vertices, vertex_id + vertices, vertex_id, vertex_id + 2 * vertices);
//            else
//          sprintf(buffer, " Vertex %d %d",
//              vertex_id, vertex_id + 2 * vertices);
//          }
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//      }
//    
//      curve_id++;
//      sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id - 4, vertex_id - 1, virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      curve_id++;
//      if(k == 1 || k == 3 || ellipse -> property == 0)
//        sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id, vertex_id - 3, virtual_curve_id);
//      else
//        sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id - 1, vertex_id - 3, virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      if(k == 1){
//        if(ellipse -> interface == 0){
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - 1 - curves, curve_id - 1, curve_id - curves, curve_id);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//        else{
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - 1 - curves, curve_id - 1 + curves, curve_id - 1, curve_id - 1 + 2 * curves);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//    
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - curves, curve_id + curves, curve_id, curve_id + 2 * curves);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//      }
//      if(k == 2){
//        sprintf(buffer, " Curve %d %d", curve_id - 1 + curves, curve_id - 1);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//        sprintf(buffer, " Curve %d %d", curve_id + curves, curve_id);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//    
//    #ifdef HUHU
//      vertex_id++;
//      sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, 
//          virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      vertex_id++;
//      if(k == 0){
//        matrix_crack_ver_id = vertex_id - 1;
//        interface_crack_ver_id = 0;
//    
//        if(ellipse -> property != 0)
//          sprintf(buffer, "vertex %d fixed vertex %d\n", vertex_id, vertex_id - 1);
//        else
//          sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, 
//              virtual_curve_id);
//      }
//      else{
//        if(k == 1){
//          matrix_crack_ver_id = vertex_id - 1;
//    
//          if(ellipse -> interface == 0)
//            sprintf(buffer, " Vertex %d %d", vertex_id - 1 - vertices, vertex_id - 1);
//          else
//            sprintf(buffer, " Vertex %d %d Vertex %d %d", 
//                vertex_id - 1 - vertices, vertex_id - 1 + vertices, vertex_id - 1, vertex_id - 1 + 2 * vertices);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//    
//          /* important: vertex fixed to vertex cannot be used in t3d2merlin;
//             parent vertex must be used instead !!! */
//          if(ellipse -> interface == 0)
//            sprintf(buffer, " Vertex %d %d", vertex_id - vertices, vertex_id - 1);
//          else
//            sprintf(buffer, " Vertex %d %d Vertex %d %d", 
//                vertex_id - 1 - vertices, vertex_id + vertices, vertex_id, vertex_id + 2 * vertices);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//    
//          sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d coincide vertex %d\n", vertex_id, 
//              crack -> center.x, crack -> center.y, 0.0, virtual_curve_id, vertex_id - 1);
//        }
//        else{
//          if(k == 3)interface_crack_ver_id = vertex_id - 1;
//    
//          sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, 
//              virtual_curve_id);
//        }
//      }
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      curve_id++;
//      sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id - 4, vertex_id - 1, virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      curve_id++;
//      sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id, vertex_id - 3, virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      if(k == 1){
//        if(ellipse -> interface == 0){
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - 1 - curves, curve_id - 1, curve_id - curves, curve_id);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//        else{
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - 1 - curves, curve_id - 1 + curves, curve_id - 1, curve_id - 1+ 2 * curves);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//    
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - curves, curve_id + curves, curve_id, curve_id + 2 * curves);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//      }
//      if(k == 2){
//        sprintf(buffer, " Curve %d %d", curve_id - 1 + curves, curve_id - 1);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//        sprintf(buffer, " Curve %d %d", curve_id + curves, curve_id);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//    #endif
//    
//      if(k == 0){
//        patch_id++;
//        if(ellipse -> property == 0)
//          sprintf(tmp_buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d", patch_id, 
//              curve_id, curve_id - 3, curve_id - 2, curve_id - 1);
//        else{
//          /*                      sprintf(buffer, "Patch %d\nElemType %d\n", patch_id, ellipse -> property); */
//          sprintf(buffer, "Patch %d\nFaceType %d\n", patch_id, ellipse -> property);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//    
//          if(imsz == 0.0){
//            default_msz = YES;
//            sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d size def property %d\n\n", patch_id, 
//                curve_id, curve_id - 3, curve_id - 2, curve_id - 1, ellipse -> property);
//          }
//          else
//            sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d size %e property %d\n\n", patch_id, 
//                curve_id, curve_id - 3, curve_id - 2, curve_id - 1, imsz, ellipse -> property);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//        }
//    
//        if(ellipse -> interface == 0){
//          if(ellipse -> property == prop || ellipse -> property == 0){
//            ellipse -> patch = patch_id;
//            break;
//          }
//        }
//        else{
//          if(ellipse -> bc != 0){
//            bc = YES;
//            sprintf(ellipse -> bc_ellipse_spec, "#Curve %d %d %d %d\n", 
//                curve_id - 3, curve_id - 2, curve_id - 1, curve_id);
//            sprintf(ellipse -> bc_matrix_spec, "#Curve %d %d %d %d\n", 
//                curve_id - 3 + curves, curve_id - 2 + curves, curve_id - 1 + curves, curve_id + curves);
//          }
//        }
//      }
//      if(k == 1){
//        patch_id++;
//        sprintf(tmp_buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d", patch_id, 
//            curve_id, curve_id - 3, curve_id - 2, curve_id - 1);
//        ellipse -> patch = patch_id;
//      }
//    
//      if(k == 1 || k == 2)write_record(ctrl_file, ctrl_file_name, "\n");
//    
//      if(k == 1 && ellipse -> interface == 0)break;
//          }
//        }
//        else{
//    
//          /* setup the points that the crack is always in the first quadrant */
//          for(i = 0; i < 4; i++){
//      g_point[i] = &(glob_point[(quadrant - 1 + i + 4) % 4]);
//      g_pnt[i] = &(glob_pnt[(quadrant - 1 + i + 4) % 4]);
//          }
//          
//          vertices = 6;
//          curves = 5 + 1;   /* one is virtual */
//    
//          for(k = 0; k < 4; k++){
//    
//      /* k = 0 inclusion boundary
//         k = 1 matrix boundary
//         k = 2 inclusion interface
//         k = 3 matrix interface */
//    
//      if(k == 1){
//        sprintf(buffer, "MasterSlave");
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//      if(k == 2){
//        sprintf(buffer, "Interface %d", ellipse -> interface);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//    
//      for(i = 0; i < 4; i++){
//        vertex_id++;
//        if(bmsz == 0.0){
//          default_msz = YES;
//          sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, g_point[i] -> x, g_point[i] -> y, 0.0);
//        }
//        else
//          sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, g_point[i] -> x, g_point[i] -> y, 0.0, bmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//        if(k == 1){
//          if(ellipse -> interface == 0)
//            sprintf(buffer, " Vertex %d %d", vertex_id - vertices, vertex_id);
//          else
//            sprintf(buffer, " Vertex %d %d Vertex %d %d", 
//                vertex_id - vertices, vertex_id + vertices, vertex_id, vertex_id + 2 * vertices);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//      }
//                  
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      for(i = 0; i < 4; i++){
//        ver_id1 = vertex_id - 3 + i;
//        ver_id2 = ver_id1 + 1;
//        if(i == 3)ver_id2 -= 4;
//                      
//        curve_id++;
//        if(i == 0){
//                          
//          /* make virtual curve over the first quadrant */
//          sprintf(buffer, "curve %d virtual order 3 vertex %d %d\n", curve_id, ver_id1, ver_id2);
//          virtual_curve_id = curve_id;
//        }
//        else{
//          if(k == 1){
//            if(ellipse -> interface == 0)
//          sprintf(buffer, " Curve %d %d", 
//              curve_id - curves, curve_id);
//            else
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - curves, curve_id + curves, curve_id, curve_id + 2 * curves);
//            write_record(ctrl_file, ctrl_file_name, buffer);
//          }
//          if(k == 2){
//            sprintf(buffer, " Curve %d %d", curve_id + curves, curve_id);
//            write_record(ctrl_file, ctrl_file_name, buffer);
//          }
//                          
//          sprintf(buffer, "curve %d order 3 vertex %d %d\n", curve_id, ver_id1, ver_id2);
//        }
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//                      
//        if(bmsz == 0.0){
//          default_msz = YES;
//          sprintf(buffer, "polygon 1 xyz %e %e %e weight %e\n", g_pnt[i] -> x, g_pnt[i] -> y, 0.0, weight_q);
//        }
//        else
//          sprintf(buffer, "polygon 1 xyz %e %e %e weight %e size %e\n", g_pnt[i] -> x, g_pnt[i] -> y, 0.0, weight_q, bmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//      }
//                  
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//          
//      /* split virtual curve */           
//      vertex_id++;
//      sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, 
//          virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      if(k == 0 || k == 1){
//        matrix_crack_ver_id = vertex_id;
//        interface_crack_ver_id = 0;
//      }
//      if(k == 3)interface_crack_ver_id = vertex_id;
//    
//      if(k == 1){
//        if(ellipse -> interface == 0)
//          sprintf(buffer, " Vertex %d %d", vertex_id - vertices, vertex_id);
//        else
//          sprintf(buffer, " Vertex %d %d Vertex %d %d", 
//              vertex_id - vertices, vertex_id + vertices, vertex_id, vertex_id + 2 * vertices);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//    
//      vertex_id++;
//      if(k == 1 || k == 3 || ellipse -> property == 0){
//        sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, 
//            virtual_curve_id);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//        if(k == 1){
//          if(ellipse -> interface == 0)
//            sprintf(buffer, " Vertex %d %d", vertex_id - 1 - vertices, vertex_id);
//          else{
//            if(ellipse -> property == 0)
//          sprintf(buffer, " Vertex %d %d Vertex %d %d",
//              vertex_id - vertices, vertex_id + vertices, vertex_id, vertex_id + 2 * vertices);
//            else
//          sprintf(buffer, " Vertex %d %d",
//              vertex_id, vertex_id + 2 * vertices);
//          }
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//      }
//    
//      curve_id++;
//      sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id - 5, vertex_id - 1, virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      curve_id++;
//      if(k == 1 || k == 3 || ellipse -> property == 0)
//        sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id, vertex_id - 4, virtual_curve_id);
//      else
//        sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id - 1, vertex_id - 4, virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      if(k == 1){
//        if(ellipse -> interface == 0){
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - 1 - curves, curve_id - 1, curve_id - curves, curve_id);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//        else{
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - 1 - curves, curve_id - 1 + curves, curve_id - 1, curve_id - 1 + 2 * curves);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//    
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - curves, curve_id + curves, curve_id, curve_id + 2 * curves);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//      }
//      if(k == 2){
//        sprintf(buffer, " Curve %d %d", curve_id - 1 + curves, curve_id - 1);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//        sprintf(buffer, " Curve %d %d", curve_id + curves, curve_id);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//    
//    #ifdef HUHU
//      vertex_id++;
//      sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, 
//          virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      vertex_id++;
//      if(k == 0){
//        matrix_crack_ver_id = vertex_id - 1;
//        interface_crack_ver_id = 0;
//    
//        if(ellipse -> property != 0)
//          sprintf(buffer, "vertex %d fixed vertex %d\n", vertex_id, vertex_id - 1);
//        else
//          sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, 
//              virtual_curve_id);
//      }
//      else{
//        if(k == 1){
//          matrix_crack_ver_id = vertex_id - 1;
//    
//          if(ellipse -> interface == 0)
//            sprintf(buffer, " Vertex %d %d", vertex_id - 1 - vertices, vertex_id - 1);
//          else
//            sprintf(buffer, " Vertex %d %d Vertex %d %d", 
//                vertex_id - 1 - vertices, vertex_id - 1 + vertices, vertex_id - 1, vertex_id - 1 + 2 * vertices);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//    
//          /* important: vertex fixed to vertex cannot be used in t3d2merlin;
//             parent vertex must be used instead !!! */
//          if(ellipse -> interface == 0)
//            sprintf(buffer, " Vertex %d %d", vertex_id - vertices, vertex_id - 1);
//          else
//            sprintf(buffer, " Vertex %d %d Vertex %d %d", 
//                vertex_id - 1 - vertices, vertex_id + vertices, vertex_id, vertex_id + 2 * vertices);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//    
//          sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d coincide vertex %d\n", vertex_id, 
//              crack -> center.x, crack -> center.y, 0.0, virtual_curve_id, vertex_id - 1);
//        }
//        else{
//          if(k == 3)interface_crack_ver_id = vertex_id - 1;
//                      
//          sprintf(buffer, "vertex %d xyz %e %e %e fixed curve %d\n", vertex_id, 
//              crack -> center.x, crack -> center.y, 0.0, virtual_curve_id);
//        }
//      }
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      curve_id++;
//      sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id - 5, vertex_id - 1, virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      curve_id++;
//      sprintf(buffer, "curve %d vertex %d %d fixed curve %d\n", curve_id, vertex_id, vertex_id - 4, virtual_curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      if(k == 1){
//        if(ellipse -> interface == 0){
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - 1 - curves, curve_id - 1, curve_id - curves, curve_id);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//        else{
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - 1 - curves, curve_id - 1 + curves, curve_id - 1, curve_id - 1+ 2 * curves);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - curves, curve_id + curves, curve_id, curve_id + 2 * curves);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//      }
//      if(k == 2){
//        sprintf(buffer, " Curve %d %d", curve_id - 1 + curves, curve_id - 1);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//        sprintf(buffer, " Curve %d %d", curve_id + curves, curve_id);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//    #endif
//    
//      if(k == 0){
//        patch_id++;
//        if(ellipse -> property == 0)
//          sprintf(tmp_buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d %d", patch_id, 
//              curve_id, curve_id - 4, curve_id - 3, curve_id - 2, curve_id - 1);
//        else{
//          /*                      sprintf(buffer, "Patch %d\nElemType %d\n", patch_id, ellipse -> property); */
//          sprintf(buffer, "Patch %d\nFaceType %d\n", patch_id, ellipse -> property);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//                          
//          if(imsz == 0)
//            sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d %d size def property %d\n\n", patch_id, 
//                curve_id, curve_id - 4, curve_id - 3, curve_id - 2, curve_id - 1, ellipse -> property);
//          else
//            sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d %d size %e property %d\n\n", patch_id, 
//                curve_id, curve_id - 4, curve_id - 3, curve_id - 2, curve_id - 1, imsz, ellipse -> property);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//        }
//    
//        if(ellipse -> interface == 0){
//          if(ellipse -> property == prop || ellipse -> property == 0){
//            ellipse -> patch = patch_id;
//            break;
//          }
//        }
//        else{
//          if(ellipse -> bc != 0){
//            bc = YES;
//            sprintf(ellipse -> bc_ellipse_spec, "#Curve %d %d %d %d %d\n", 
//                curve_id - 4, curve_id - 3, curve_id - 2, curve_id - 1, curve_id);
//            sprintf(ellipse -> bc_matrix_spec, "#Curve %d %d %d %d %d\n", 
//                curve_id - 4 + curves, curve_id - 3 + curves, curve_id - 2 + curves, curve_id - 1 + curves, curve_id + curves);
//          }
//        }
//      }
//      if(k == 1){
//        patch_id++;
//        sprintf(tmp_buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d %d", patch_id, 
//            curve_id, curve_id - 4, curve_id - 3, curve_id - 2, curve_id - 1);
//        ellipse -> patch = patch_id;
//      }
//    
//      if(k == 1 || k == 2)write_record(ctrl_file, ctrl_file_name, "\n");
//    
//      if(k == 1 && ellipse -> interface == 0)break;
//          }
//        }
//    
//        /* make crack */
//        if(ellipse -> interface == 0){
//    
//          /* make new matrix crack vertices because the original ones have got a weight != 1.0, 
//       while the interface vertices (to be created later because interface_crack_ver_id == 0) will have weight == 1.0;
//       this may cause non-compatible meshes on master/slave curves */
//          vertex_id++;
//          sprintf(buffer, "vertex %d fixed vertex %d weight 1.0\n", vertex_id, matrix_crack_ver_id);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//          vertex_id++;
//          if(ellipse -> property != prop)
//      sprintf(buffer, "vertex %d fixed vertex %d weight 1.0\n", vertex_id, matrix_crack_ver_id + 1);
//          else
//      sprintf(buffer, "vertex %d fixed vertex %d weight 1.0\n", vertex_id, matrix_crack_ver_id);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//          orig_matrix_crack_ver_id = matrix_crack_ver_id;
//          matrix_crack_ver_id = vertex_id - 1;
//        }
//    
//        vertex_id++;
//        if(cmsz == 0.0){
//          default_msz = YES;
//          sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0);
//        }
//        else
//          sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0, cmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//        tip_crack_ver_id = vertex_id;
//    
//        write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//        /* note: cracksurface, interface and master/slave curves must be of the same orientation;
//           to enable contourpath, the curves must go from tip to center !!! */
//        curve_id++;
//        if(cmsz == 0.0)
//          sprintf(buffer, "curve %d vertex %d %d\n", curve_id, vertex_id, matrix_crack_ver_id);
//        else
//          sprintf(buffer, "curve %d vertex %d %d size %e\n", curve_id, vertex_id, matrix_crack_ver_id, cmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//        crack_cur_id = curve_id;
//    
//        curve_id++;
//        if(cmsz == 0.0)
//          sprintf(buffer, "curve %d vertex %d %d coincide curve %d\n", curve_id, vertex_id, matrix_crack_ver_id + 1, curve_id - 1);
//        else
//          sprintf(buffer, "curve %d vertex %d %d coincide curve %d size %e\n", curve_id, vertex_id, matrix_crack_ver_id + 1,
//            curve_id - 1, cmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//        write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//        /* make crack interface */
//        if(crack -> interface != 0){
//          if(ellipse -> interface == 0){
//      vertex_id++;
//      if(cmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, crack -> center.x, crack -> center.y, 0.0);
//      }
//      else
//        sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      vertex_id++;
//      if(cmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, crack -> center.x, crack -> center.y, 0.0);
//      }
//      else
//        sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      interface_crack_ver_id = vertex_id - 1;
//          }
//              
//          vertex_id++;
//          if(cmsz == 0.0){
//      default_msz = YES;
//      sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0);
//          }
//          else
//      sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0, cmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//              
//          vertex_id++;
//          if(cmsz == 0.0){
//      default_msz = YES;
//      sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0);
//          }
//          else
//      sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0, cmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//              
//          write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//          /* note: cracksurface, interface and master/slave curves must be of the same orientation;
//       to enable contourpath, the curves must go from tip to center !!! */
//          curve_id++;
//          if(cmsz == 0.0)
//      sprintf(buffer, "curve %d vertex %d %d\n", curve_id, vertex_id - 1, interface_crack_ver_id);
//          else
//      sprintf(buffer, "curve %d vertex %d %d size %e\n", curve_id, vertex_id - 1, interface_crack_ver_id, cmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//          curve_id++;
//          if(cmsz == 0.0)
//      sprintf(buffer, "curve %d vertex %d %d\n", curve_id, vertex_id, interface_crack_ver_id + 1);
//          else
//      sprintf(buffer, "curve %d vertex %d %d size %e\n", curve_id, vertex_id, interface_crack_ver_id + 1, cmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//          write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//          if(ellipse -> interface == 0){
//    
//      /* important: vertex fixed to vertex cannot be used in t3d2merlin;
//         parent vertex must be used instead !!! */
//      if(ellipse -> property != 0){
//        if(ellipse -> property != prop)
//          sprintf(buffer, "MasterSlave Vertex %d %d Vertex %d %d Vertex %d %d Vertex %d %d Curve %d %d Curve %d %d\n",
//              tip_crack_ver_id, vertex_id, tip_crack_ver_id, vertex_id - 1, 
//              orig_matrix_crack_ver_id - vertices, interface_crack_ver_id, 
//              orig_matrix_crack_ver_id - vertices, interface_crack_ver_id + 1, 
//              curve_id - 3, curve_id - 1, curve_id - 2, curve_id);
//        else
//          sprintf(buffer, "MasterSlave Vertex %d %d Vertex %d %d Vertex %d %d Vertex %d %d Curve %d %d Curve %d %d\n", 
//              tip_crack_ver_id, vertex_id, tip_crack_ver_id, vertex_id - 1, 
//              orig_matrix_crack_ver_id, interface_crack_ver_id, orig_matrix_crack_ver_id, interface_crack_ver_id + 1, 
//              curve_id - 3, curve_id - 1, curve_id - 2, curve_id);
//      }
//      else
//        sprintf(buffer, "MasterSlave Vertex %d %d Vertex %d %d Vertex %d %d Vertex %d %d Curve %d %d Curve %d %d\n", 
//            tip_crack_ver_id, vertex_id, tip_crack_ver_id, vertex_id - 1, 
//            orig_matrix_crack_ver_id, interface_crack_ver_id, orig_matrix_crack_ver_id + 1, interface_crack_ver_id + 1, 
//            curve_id - 3, curve_id - 1, curve_id - 2, curve_id);
//          }
//          else
//      sprintf(buffer, "MasterSlave Vertex %d %d Vertex %d %d Curve %d %d Curve %d %d\n", 
//          tip_crack_ver_id, vertex_id, tip_crack_ver_id, vertex_id - 1, 
//          curve_id - 3, curve_id - 1, curve_id - 2, curve_id);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//          sprintf(buffer, "Interface %d Curve %d %d\n", crack -> interface, curve_id, curve_id - 1);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//    
//        crack_id++;
//        sprintf(buffer, "CrackSurface %d Curve %d %d\n", crack_id, crack_cur_id, crack_cur_id + 1);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//        sprintf(buffer, "Vertex %d\nCrackFronts %d\n", tip_crack_ver_id, crack_id);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//          
//        if(crack -> bc != 0){
//          bc = YES;
//          sprintf(crack -> bc_spec, "#Curve %d %d\n", crack_cur_id, crack_cur_id + 1);
//        }
//    
//        if(crack -> radius != 0.0){
//          sprintf(buffer, "Contourpath %d Radius %e\n", crack_id, crack -> radius);
//          write_record(ctrl_file, ctrl_file_name, buffer);
//        }
//    
//        if(tmsz != 0.0){
//          vertex_id++;
//          sprintf(buffer, "vertex %d virtual fixed vertex %d size %e\n", vertex_id, tip_crack_ver_id, tmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//          write_record(t3d_in_file, t3d_in_file_name, "\n");
//        }
//    
//        if(ellipse -> property == 0 || ellipse -> interface != 0 || ellipse -> property != prop){
//          sprintf(buffer, "%s %d %d hole\n\n", tmp_buffer, -crack_cur_id, crack_cur_id + 1);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//        }
//        else
//          crack -> curve = crack_cur_id;
//      }
//      else{
//        strcat(buffer, "\n\n");
//    
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//    
//        vertices = 4;
//        curves = 4;
//    
//        for(k = 0; k < 4; k++){
//    
//          /* k = 0 inclusion boundary
//       k = 1 matrix boundary
//       k = 2 inclusion interface
//       k = 3 matrix interface */
//    
//          if(k == 1){
//      sprintf(buffer, "MasterSlave");
//      write_record(ctrl_file, ctrl_file_name, buffer);
//          }
//          if(k == 2){
//      sprintf(buffer, "Interface %d", ellipse -> interface);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//          }
//    
//          for(i = 0; i < 4; i++){
//      vertex_id++;
//      if(bmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, glob_point[i].x, glob_point[i].y, 0.0);
//      }
//      else
//        sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, glob_point[i].x, glob_point[i].y, 0.0, bmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      if(k == 1){
//        if(ellipse -> interface == 0)
//          sprintf(buffer, " Vertex %d %d", vertex_id - vertices, vertex_id);
//        else
//          sprintf(buffer, " Vertex %d %d Vertex %d %d", 
//              vertex_id - vertices, vertex_id + vertices, vertex_id, vertex_id + 2 * vertices);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//          }
//          
//          write_record(t3d_in_file, t3d_in_file_name, "\n");
//              
//          for(i = 0; i < 4; i++){
//      ver_id1 = vertex_id - 3 + i;
//      ver_id2 = ver_id1 + 1;
//      if(i == 3)ver_id2 -= 4;
//                  
//      curve_id++;
//      sprintf(buffer, "curve %d order 3 vertex %d %d\n", curve_id, ver_id1, ver_id2);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//                  
//      if(bmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "polygon 1 xyz %e %e %e weight %e\n", glob_pnt[i].x, glob_pnt[i].y, 0.0, weight_q);
//      }
//      else
//        sprintf(buffer, "polygon 1 xyz %e %e %e weight %e size %e\n", glob_pnt[i].x, glob_pnt[i].y, 0.0, weight_q, bmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      if(k == 1){
//        if(ellipse -> interface == 0)
//          sprintf(buffer, " Curve %d %d", curve_id - curves, curve_id);
//        else
//          sprintf(buffer, " Curve %d %d Curve %d %d", 
//              curve_id - curves, curve_id + curves, curve_id, curve_id + 2 * curves);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//      if(k == 2){
//        sprintf(buffer, " Curve %d %d", curve_id + curves, curve_id);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//          }
//          
//          write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//          if(k == 0){
//      patch_id++;
//      if(ellipse -> property == 0)
//        sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d hole\n\n", patch_id, 
//            curve_id - 3, curve_id - 2, curve_id - 1, curve_id);
//      else{
//        /*                    sprintf(buffer, "Patch %d\nElemType %d\n", patch_id, ellipse -> property); */
//        sprintf(buffer, "Patch %d\nFaceType %d\n", patch_id, ellipse -> property);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//    
//        if(imsz == 0.0){
//          default_msz = YES;
//          sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d size def property %d\n\n", patch_id, 
//              curve_id - 3, curve_id - 2, curve_id - 1, curve_id, ellipse -> property);
//        }
//        else
//          sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d size %e property %d\n\n", patch_id, 
//              curve_id - 3, curve_id - 2, curve_id - 1, curve_id, imsz, ellipse -> property);
//      }
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      if(ellipse -> interface == 0){
//        if(ellipse -> property == prop || ellipse -> property == 0){
//          ellipse -> patch = patch_id;
//          break;
//        }
//      }
//      else{
//        if(ellipse -> bc != 0){
//          bc = YES;
//          sprintf(ellipse -> bc_ellipse_spec, "#Curve %d %d %d %d\n", 
//              curve_id - 3, curve_id - 2, curve_id - 1, curve_id);
//          sprintf(ellipse -> bc_matrix_spec, "#Curve %d %d %d %d\n", 
//              curve_id - 3 + curves, curve_id - 2 + curves, curve_id - 1 + curves, curve_id + curves);
//        }
//      }
//          }
//          if(k == 1){
//      patch_id++;
//      sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d hole\n\n", patch_id, 
//          curve_id - 3, curve_id - 2, curve_id - 1, curve_id);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      ellipse -> patch = patch_id;
//          }
//    
//          if(k == 1 || k == 2)write_record(ctrl_file, ctrl_file_name, "\n");
//    
//          if(k == 1 && ellipse -> interface == 0)break;
//        }
//      }
//    }
//    
//    
//    
//    static void
//    write_crack(crack_rec *crack)
//    {
//      point_rec mid_pnt;
//      double cmsz, tmsz;
//      int i, tip_crack_ver_id1, tip_crack_ver_id2, mid_crack_ver_id1, mid_crack_ver_id2, crack_cur_id1, crack_cur_id2;
//      char buffer[1024];
//    
//      if(crack -> id < 0)return;
//      if(crack -> ellipse != NULL)return;
//    
//      if(crack -> interface == 0)
//        sprintf(buffer, "\n# crack no %d\n\n", crack -> id);
//      else
//        sprintf(buffer, "\n# crack no %d (interface %d)\n\n", crack -> id, crack -> interface);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//    
//      cmsz = crack -> crack_msz;
//      tmsz = crack -> tip_msz;
//    
//      /* note: separate crack is modelled as two cracks (oriented from tip to mid_pnt) */
//    
//      aver_point(mid_pnt, crack -> center, crack -> tip);
//    
//      /* make crack */
//      vertex_id++;
//      if(cmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, crack -> center.x, crack -> center.y, 0.0);
//      }
//      else
//        sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      tip_crack_ver_id1 = vertex_id;
//      
//      vertex_id++;
//      if(cmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0);
//      }
//      else
//        sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      tip_crack_ver_id2 = vertex_id;
//      
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      vertex_id++;
//      if(cmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, mid_pnt.x, mid_pnt.y, 0.0);
//      }
//      else
//        sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, mid_pnt.x, mid_pnt.y, 0.0, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      mid_crack_ver_id1 = vertex_id;
//      
//      vertex_id++;
//      if(cmsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "vertex %d xyz %e %e %e coincide vertex %d\n", vertex_id, mid_pnt.x, mid_pnt.y, 0.0, vertex_id - 1);
//      }
//      else
//        sprintf(buffer, "vertex %d xyz %e %e %e coincide vertex %d size %e\n", vertex_id, mid_pnt.x, mid_pnt.y, 0.0, 
//          vertex_id - 1, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      mid_crack_ver_id2 = vertex_id;
//      
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      /* make crack interface */
//      if(crack -> interface != 0){
//        for(i = 0; i < 2; i++){
//          vertex_id++;
//          if(cmsz == 0.0){
//      default_msz = YES;
//      sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, crack -> center.x, crack -> center.y, 0.0);
//          }
//          else
//      sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, crack -> center.x, crack -> center.y, 0.0, cmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//        }
//          
//        for(i = 0; i < 2; i++){
//          vertex_id++;
//          if(cmsz == 0.0){
//      default_msz = YES;
//      sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0);
//          }
//          else
//      sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, crack -> tip.x, crack -> tip.y, 0.0, cmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//        }
//      
//        write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//        for(i = 0; i < 2; i++){
//          vertex_id++;
//          if(cmsz == 0.0){
//      default_msz = YES;
//      sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, mid_pnt.x, mid_pnt.y, 0.0);
//          }
//          else
//      sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, mid_pnt.x, mid_pnt.y, 0.0, cmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//        }
//      
//        write_record(t3d_in_file, t3d_in_file_name, "\n");
//      }
//    
//      /* make crack */
//      /* note: cracksurface, interface and master/slave curves must be of the same orientation (from tip to mid_pnt) */
//      curve_id++;
//      if(cmsz == 0.0)
//        sprintf(buffer, "curve %d vertex %d %d\n", curve_id, tip_crack_ver_id1, mid_crack_ver_id1);
//      else
//        sprintf(buffer, "curve %d vertex %d %d size %e\n", curve_id, tip_crack_ver_id1, mid_crack_ver_id1, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      
//      crack_cur_id1 = curve_id;
//    
//      curve_id++;
//      if(cmsz == 0.0)
//        sprintf(buffer, "curve %d vertex %d %d coincide curve %d\n", curve_id, tip_crack_ver_id1, mid_crack_ver_id2, curve_id - 1);
//      else
//        sprintf(buffer, "curve %d vertex %d %d coincide curve %d size %e\n", curve_id, tip_crack_ver_id1, mid_crack_ver_id2, 
//          curve_id - 1, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      
//      curve_id++;
//      if(cmsz == 0.0)
//        sprintf(buffer, "curve %d vertex %d %d\n", curve_id, tip_crack_ver_id2, mid_crack_ver_id1);
//      else
//        sprintf(buffer, "curve %d vertex %d %d size %e\n", curve_id, tip_crack_ver_id2, mid_crack_ver_id1, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      crack_cur_id2 = curve_id;
//      
//      curve_id++;
//      if(cmsz == 0.0)
//        sprintf(buffer, "curve %d vertex %d %d coincide curve %d\n", curve_id, tip_crack_ver_id2, mid_crack_ver_id2, curve_id - 1);
//      else
//        sprintf(buffer, "curve %d vertex %d %d coincide curve %d size %e\n", curve_id, tip_crack_ver_id2, mid_crack_ver_id2, 
//          curve_id - 1, cmsz);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//      
//      /* make crack interface */
//      /* note: cracksurface, interface and master/slave curves must be of the same orientation (from tip to mid_pnt) */
//      if(crack -> interface != 0){
//        for(i = 0; i < 2; i++){
//          curve_id++;
//          if(cmsz == 0.0)
//      sprintf(buffer, "curve %d vertex %d %d\n", curve_id, vertex_id - 5 + i, vertex_id - 1 + i);
//          else
//      sprintf(buffer, "curve %d vertex %d %d size %e\n", curve_id, vertex_id - 5 + i, vertex_id - 1 + i, cmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//        }
//        for(i = 0; i < 2; i++){
//          curve_id++;
//          if(cmsz == 0.0)
//      sprintf(buffer, "curve %d vertex %d %d\n", curve_id, vertex_id - 3 + i, vertex_id - 1 + i);
//          else
//      sprintf(buffer, "curve %d vertex %d %d size %e\n", curve_id, vertex_id - 3 + i, vertex_id - 1 + i, cmsz);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//        }
//      
//        sprintf(buffer, "MasterSlave Vertex %d %d Vertex %d %d Vertex %d %d Vertex %d %d Vertex %d %d Vertex %d %d", 
//          tip_crack_ver_id1, vertex_id - 5, tip_crack_ver_id1, vertex_id - 4, 
//          tip_crack_ver_id2, vertex_id - 3, tip_crack_ver_id2, vertex_id - 2, 
//          mid_crack_ver_id1, vertex_id - 1, mid_crack_ver_id2, vertex_id);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//        sprintf(buffer, " Curve %d %d Curve %d %d Curve %d %d Curve %d %d\n", 
//          curve_id - 7, curve_id - 3, curve_id - 6, curve_id - 2,
//          curve_id - 5, curve_id - 1, curve_id - 4, curve_id);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//        sprintf(buffer, "Interface %d Curve %d %d Curve %d %d\n", crack -> interface, curve_id - 3, curve_id - 2,
//          curve_id, curve_id - 1);
//        write_record(ctrl_file, ctrl_file_name, buffer);
//      }
//      
//      crack_id++;
//      sprintf(buffer, "CrackSurface %d Curve %d %d\n", crack_id, crack_cur_id1, crack_cur_id1 + 1);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//      sprintf(buffer, "Vertex %d\nCrackFronts %d\n", tip_crack_ver_id1, crack_id);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//      sprintf(buffer, "Contourpath %d Radius %e\n", crack_id, crack -> radius);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//      
//      crack_id++;
//      sprintf(buffer, "CrackSurface %d Curve %d %d\n", crack_id, crack_cur_id2, crack_cur_id2 + 1);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//      sprintf(buffer, "Vertex %d\nCrackFronts %d\n", tip_crack_ver_id2, crack_id);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//      sprintf(buffer, "Contourpath %d Radius %e\n", crack_id, crack -> radius);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//          
//      if(crack -> bc != 0){
//        bc = YES;
//        sprintf(crack -> bc_spec, "#Curve %d %d %d %d\n", crack_cur_id1, crack_cur_id1 + 1, crack_cur_id2, crack_cur_id2 + 1);
//      }
//      
//      if(tmsz != 0.0){
//        vertex_id++;
//        sprintf(buffer, "vertex %d virtual fixed vertex %d size %e\n", vertex_id, tip_crack_ver_id1, tmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//        vertex_id++;
//        sprintf(buffer, "vertex %d virtual fixed vertex %d size %e\n", vertex_id, tip_crack_ver_id2, tmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//          
//        write_record(t3d_in_file, t3d_in_file_name, "\n");
//      }
//    
//      crack -> curve = crack_cur_id1;
//    }
//    
//    
//    
//    static void
//    write_cylinder(cylinder_rec *cylinder, ellipse_rec *ellipse_array, crack_rec *crack_array, int ellipses, int cracks)
//    {
//      point_rec *point[4];
//      ellipse_rec *ellipse = NULL;
//      crack_rec *crack = NULL;
//      double bmsz, imsz;
//      int i, ver_id1, ver_id2, ellipse_count, crack_count, num;
//      char buffer[1024];
//    
//      bmsz = cylinder -> boundary_msz;
//      imsz = cylinder -> internal_msz;
//    
//      sprintf(buffer, "\n# matrix (material %d)\n\n", cylinder -> property);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//    
//      point[0] = &(cylinder -> lower_left);
//      point[1] = &(cylinder -> lower_right);
//      point[2] = &(cylinder -> upper_right);
//      point[3] = &(cylinder -> upper_left);
//    
//      for(i = 0; i < 4; i++){
//        vertex_id++;
//        if(bmsz == 0.0){
//          default_msz = YES;
//          sprintf(buffer, "vertex %d xyz %e %e %e\n", vertex_id, point[i] -> x, point[i] -> y, 0.0);
//        }
//        else
//          sprintf(buffer, "vertex %d xyz %e %e %e size %e\n", vertex_id, point[i] -> x, point[i] -> y, 0.0, bmsz);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//      }
//    
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      for(i = 0; i < 4; i++){
//        ver_id1 = vertex_id - 3 + i;
//        ver_id2 = ver_id1 + 1;
//        if(i == 3)ver_id2 -= 4;
//    
//        curve_id++;
//        sprintf(buffer, "curve %d vertex %d %d\n", curve_id, ver_id1, ver_id2);
//        write_record(t3d_in_file, t3d_in_file_name, buffer);
//      }
//    
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      patch_id++;
//      if(imsz == 0.0){
//        default_msz = YES;
//        sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d size def property %d", patch_id, 
//          curve_id - 3, curve_id - 2, curve_id - 1, curve_id, cylinder -> property);
//      }
//      else
//        sprintf(buffer, "patch %d normal 0.0 0.0 1.0 boundary curve %d %d %d %d size %e property %d", patch_id, 
//          curve_id - 3, curve_id - 2, curve_id - 1, curve_id, imsz, cylinder -> property);
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//    
//      /*  sprintf(buffer, "Patch %d\nElemType %d\n", patch_id, cylinder -> property); */
//      sprintf(buffer, "Patch %d\nFaceType %d\n", patch_id, cylinder -> property);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//    
//      ellipse_count = 0;
//      for(i = 0; i < ellipses; i++){
//        ellipse = &(ellipse_array[i]);
//        if(ellipse -> patch == 0)continue;
//        ellipse_count++;
//      }
//    
//      if(ellipse_count != 0){
//        write_record(t3d_in_file, t3d_in_file_name, " \\\nsubpatch");
//    
//        num = 0;
//        for(i = 0; i < ellipses; i++){
//          ellipse = &(ellipse_array[i]);
//          if(ellipse -> patch == 0)continue;
//          num++;
//    
//          sprintf(buffer, " %d", ellipse -> patch);
//          write_record(t3d_in_file, t3d_in_file_name, buffer);
//              
//          if(num % 30 == 0 && ellipse_count - num > 5)
//      write_record(t3d_in_file, t3d_in_file_name, " \\\n        ");
//        }     
//      }
//    
//      crack_count = 0;
//      for(i = 0; i < cracks; i++){
//        crack = &(crack_array[i]);
//        if(crack -> curve == 0)continue;
//        crack_count += 2;
//        if(crack -> ellipse == NULL)crack_count += 2;
//      }
//    
//      if(crack_count != 0){
//        write_record(t3d_in_file, t3d_in_file_name, " \\\nboundary curve");
//          
//        num = 0;
//        for(i = 0; i < cracks; i++){
//          crack = &(crack_array[i]);
//          if(crack -> curve == 0)continue;
//          num++;
//    
//          if(crack -> ellipse == NULL){
//      num++;
//      sprintf(buffer, " %d %d %d %d", -crack -> curve, crack -> curve + 1, crack -> curve + 2, -(crack -> curve + 3));
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//          }
//          else{
//      sprintf(buffer, " %d %d", crack -> curve, -(crack -> curve + 1));
//      write_record(t3d_in_file, t3d_in_file_name, buffer);
//          }
//    
//          if(num % 10 == 0 && crack_count - num > 2)
//      write_record(t3d_in_file, t3d_in_file_name, " \\\n              ");
//        }     
//      }
//    
//      write_record(t3d_in_file, t3d_in_file_name, "\n");
//    
//      sprintf(buffer, "\n# Vertex %d (bottom left) %d (bottom right) %d (top right) %d (top left)\n",
//        vertex_id - 3, vertex_id - 2, vertex_id - 1, vertex_id);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//      sprintf(buffer, "# Curve %d (bottom) %d (right) %d (top) %d (left)\n\n",
//        curve_id - 3, curve_id - 2, curve_id - 1, curve_id);
//      write_record(ctrl_file, ctrl_file_name, buffer);
//    }
//          
//          
//    
//    static logic
//    check_ellipse_inside_cylinder(ellipse_rec *ellipse, cylinder_rec *cylinder)
//    {
//      double cx, cy, max, min;
//    
//      cx = ellipse -> center.x;
//      cy = ellipse -> center.y;
//    
//      max = ellipse -> max + epsilon;
//      min = ellipse -> min + epsilon;
//    
//      if(cx - min_x <= min || max_x - cx <= min)return(NO);
//      if(cy - min_y <= min || max_y - cy <= min)return(NO);
//    
//      if(cx - min_x <= max){
//        if(check_ellipse_line_intersection(ellipse, &(cylinder -> lower_left), &(cylinder -> upper_left)) == YES)return(NO);
//      }
//      if(max_x - cx <= max){
//        if(check_ellipse_line_intersection(ellipse, &(cylinder -> lower_right), &(cylinder -> upper_right)) == YES)return(NO);
//      }
//      if(cy - min_y <= max){
//        if(check_ellipse_line_intersection(ellipse, &(cylinder -> lower_left), &(cylinder -> lower_right)) == YES)return(NO);
//      }
//      if(max_y - cy <= max){
//        if(check_ellipse_line_intersection(ellipse, &(cylinder -> upper_left), &(cylinder -> upper_right)) == YES)return(NO);
//      }
//    
//      return(YES);
//    }
//    
//    
//

bool check_ellipse_ellipse_overlap (ellipse_rec *ellipse1, ellipse_rec *ellipse2)
{
  point_rec pnt1, pnt2, p1, p2;
  double max, max1, max2;
  double dist, dist1, dist2;
  double t1, tt1, ttt1, t2, tt2, ttt2, delta_t = EPSILON_T;
  int qq1, qqq1, qq2, qqq2;
  //int q1, q2;

  if (ellipse1 -> id == ellipse2 -> id)return(NO);

  max1 = ellipse1 -> max + epsilon;
  max2 = ellipse2 -> max + epsilon;
  max = max1 + max2;

  /* intersection of circumsribed circles */
  dist = dist_point (ellipse1 -> center, ellipse2 -> center);
  if (dist > max * max)return(NO);

  /* center of ellipse2 inside ellipse1 */
  if (check_point_inside_ellipse (&(ellipse2 -> center), ellipse1) == YES)return(YES);

  /* center of ellipse1 inside ellipse2 */
  if (check_point_inside_ellipse (&(ellipse1 -> center), ellipse2) == YES)return(YES);

  /* intersection of circle circumscribed to ellipse1 with ellipse2 */
  //q2 =
  project_point_to_ellipse (&(ellipse1 -> center), ellipse2, &pnt2, &t2);
  dist1 = dist_point (ellipse1 -> center, pnt2);
  if (dist1 > max1 * max1)return(NO);

  /* closest point on ellipse2 to center of ellipse1 inside ellipse1 */
  if (check_point_inside_ellipse (&pnt2, ellipse1) == YES)return(YES);

  /* intersection of circle circumscribed to ellipse2 with ellipse1 */
  //q1 =
  project_point_to_ellipse (&(ellipse2 -> center), ellipse1, &pnt1, &t1);
  dist2 = dist_point (ellipse2 -> center, pnt1);
  if (dist2 > max2 * max2)return(NO);

  /* closest point on ellipse1 to center of ellipse2 inside ellipse2 */
  if (check_point_inside_ellipse (&pnt1, ellipse2) == YES)return(YES);

  /* the iterations should not oscillate between quadrants because the circle tests have been already done */

  // termitovo: musel jsem zmensit delta_t, jinak to nedostatecne presne iterovalo
  // printf mi ukaze jak casto se to sem dostane
  //printf ("\nmeso3d iteration\n"); fflush(stdin);

  /* find closest points starting from pnt1 */
  while (YES) {
    qqq2 = project_point_to_ellipse (&pnt1, ellipse2, &p2, &ttt2);
    qq1 = project_point_to_ellipse (&p2, ellipse1, &pnt1, &tt1);
    if (fabs (t1 - tt1) < delta_t)break;
    t1 = tt1;
  }

  /* find closest points starting from pnt2 */
  while (YES) {
    qqq1 = project_point_to_ellipse (&pnt2, ellipse1, &p1, &ttt1);
    qq2 = project_point_to_ellipse (&p1, ellipse2, &pnt2, &tt2);
    if (fabs (t2 - tt2) < delta_t)break;
    t2 = tt2;
  }

  /* there is a problem that intersections of ellipses may be in different quadrants;
     in this case averaging of parameters makes no sense but it is clear that ellipses do intersect */

  if (qq1 != qqq1 || qq2 != qqq2)return(YES);

  t1 = (tt1 + ttt1) / 2.0;
  t2 = (tt2 + ttt2) / 2.0;
  get_ellipse_point_quadrant (ellipse1, qq1, t1, &pnt1);
  get_ellipse_point_quadrant (ellipse2, qq2, t2, &pnt2);

  /* check pnt1 inside ellipse2 */
  if (check_point_inside_ellipse (&pnt1, ellipse2) == YES)return(YES);

  /* check pnt2 inside ellipse1 */
  if (check_point_inside_ellipse (&pnt2, ellipse1) == YES)return(YES);

  return(NO);
}


//    
//    
//    
//    static logic
//    check_ellipse_line_intersection(ellipse_rec *ellipse, point_rec *point1, point_rec *point2)
//    {
//      point_rec pnt1, pnt2;
//      double a0, a1, a2, a0x, a0y, a1x, a1y, a2x, a2y, D, max, min, t1, t2, dx, dy, dd;
//    
//      max = ellipse -> max + epsilon;
//      min = ellipse -> min + epsilon;
//    
//      transform_from_global_to_local(point1, &pnt1, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//      transform_from_global_to_local(point2, &pnt2, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//    
//      dx = pnt2.x - pnt1.x;
//      dy = pnt2.y - pnt1.y;
//      dd = sqrt(dx * dx + dy * dy);
//    
//      if(dd < EPSILON)error_message(GENERAL_ERROR, "Inclusion %d checked against intersection with zero length line", ellipse -> id);
//    
//      a0x = pnt1.x / max;
//      a0y = pnt1.y / min;
//      a0 = a0x * a0x + a0y * a0y - 1.0;
//    
//      a2x = dx / max;
//      a2y = dy / min;
//      a2 = a2x * a2x + a2y * a2y;
//    
//      a1x = a0x * a2x;
//      a1y = a0y * a2y;
//      a1 = a1x + a1y;
//    
//      D = a1 * a1 - a0 * a2;
//      if(D < 0.0)return(NO);
//      D = sqrt(D);
//    
//      t1 = -(a1 - D) / a2;
//      t2 = -(a1 + D) / a2;
//    
//      if(t1 >= 0.0 && t1 <= 1.0)return(YES);
//      if(t2 >= 0.0 && t2 <= 1.0)return(YES);
//    
//      return(NO);
//    }
//    
//    
//    
//    static logic
//    check_crack_cylinder_intersection(crack_rec *crack, cylinder_rec *cylinder)
//    {
//      if(check_line_line_intersection(&(crack -> center), &(crack -> tip), 
//                    &(cylinder -> lower_left), &(cylinder -> lower_right)) == YES)return(YES);
//      if(check_line_line_intersection(&(crack -> center), &(crack -> tip), 
//                    &(cylinder -> upper_left), &(cylinder -> upper_right)) == YES)return(YES);
//      if(check_line_line_intersection(&(crack -> center), &(crack -> tip), 
//                    &(cylinder -> lower_left), &(cylinder -> upper_left)) == YES)return(YES);
//      if(check_line_line_intersection(&(crack -> center), &(crack -> tip), 
//                    &(cylinder -> lower_right), &(cylinder -> upper_right)) == YES)return(YES);
//      return(NO);
//    }
//    
//    
//    
//    static logic
//    check_crack_ellipse_intersection(crack_rec *crack, ellipse_rec *ellipse)
//    {
//      if(crack -> ellipse == ellipse)return(NO);
//    
//      return(check_ellipse_line_intersection(ellipse, &(crack -> center), &(crack -> tip)));
//    }
//    
//    
//    
//    static logic
//    check_crack_crack_intersection(crack_rec *crack1, crack_rec *crack2)
//    {
//      if(crack1 -> id == crack2 -> id)return(NO);
//    
//      return(check_line_line_intersection(&(crack1 -> center), &(crack1 -> tip), &(crack2 -> center), &(crack2 -> tip)));
//    }
//    
//    
//    
//    static logic
//    check_crack_inside_ellipse(crack_rec *crack, ellipse_rec *ellipse)
//    {
//      if(check_point_inside_ellipse(&(crack -> center), ellipse) == YES)return(YES);
//      if(check_point_inside_ellipse(&(crack -> tip), ellipse) == YES)return(YES);
//    
//      return(NO);
//    }
//    
//    
//

static logic check_point_inside_ellipse (point_rec *point, ellipse_rec *ellipse)
{
  point_rec pnt;
  double x_rate, y_rate;

  transform_from_global_to_local (point, &pnt, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);

  x_rate = pnt.x / (ellipse -> max + epsilon);
  y_rate = pnt.y / (ellipse -> min + epsilon);
  if (x_rate * x_rate + y_rate * y_rate <= 1.0)return(YES);

  return(NO);
}

//    
//    
//    
//    static logic
//    check_line_line_intersection(point_rec *pnta1, point_rec *pnta2, point_rec *pntb1, point_rec *pntb2)
//    {
//      point_rec pointa1, pointa2, pointb1, pointb2;
//      double dax, day, dbx, dby, dx, dy, ta, tb, dda, ddb, dd, dist, product, det;
//    
//      dax = pnta2 -> x - pnta1 -> x;
//      day = pnta2 -> y - pnta1 -> y;
//    
//      dbx = pntb2 -> x - pntb1 -> x;
//      dby = pntb2 -> y - pntb1 -> y;
//    
//      dda = sqrt(dax * dax + day * day);
//      ddb = sqrt(dbx * dbx + dby * dby);
//    
//      if(dda < EPSILON)error_message(GENERAL_ERROR, "Line intersection with zero length line");
//      if(ddb < EPSILON)error_message(GENERAL_ERROR, "Line intersection with zero length line");
//    
//      if(epsilon != 0){
//    
//        /* shift points by epsilon on each side */
//        pointa1.x = pnta1 -> x - dax / dda * epsilon;
//        pointa1.y = pnta1 -> y - day / dda * epsilon;
//        pointa2.x = pnta2 -> x - dax / dda * epsilon;
//        pointa2.y = pnta2 -> y - day / dda * epsilon;
//    
//        pointb1.x = pntb1 -> x - dbx / ddb * epsilon;
//        pointb1.y = pntb1 -> y - dby / ddb * epsilon;
//        pointb2.x = pntb2 -> x - dbx / ddb * epsilon;
//        pointb2.y = pntb2 -> y - dby / ddb * epsilon;
//          
//        dda += 2.0 * epsilon;
//        ddb += 2.0 * epsilon;
//    
//        pnta1 = &pointa1;
//        pnta2 = &pointa2;
//    
//        pntb1 = &pointb1;
//        pntb2 = &pointb2;
//      }
//          
//      dx = pntb1 -> x - pnta1 -> x;
//      dy = pntb1 -> y - pnta1 -> y;
//      
//      det = dbx * day - dby * dax;
//    
//      /* check for parallel lines */
//      if(fabs(det) < EPSILON){
//    
//        /* check their distance   */
//        dd = dx * dx + dy * dy;
//        if(dd < epsilon * epsilon || dd < EPSILON * EPSILON)return(YES);
//    
//        product = (dax * dx + day * dy) / dda;
//        dist = fabs(1.0 - product * product / dd); 
//        if(dist > epsilon * epsilon && dist > EPSILON * EPSILON)return(NO);
//    
//        /* check if they overlap */
//        dist = dd;
//        dx = pntb2 -> x - pnta2 -> x;
//        dy = pntb2 -> y - pnta2 -> y;
//        dd = dx * dx + dy * dy;
//        if(dist < dd)dist = dd;
//        dx = pntb2 -> x - pnta1 -> x;
//        dy = pntb2 -> y - pnta1 -> y;
//        dd = dx * dx + dy * dy;
//        if(dist < dd)dist = dd;
//        dx = pntb1 -> x - pnta2 -> x;
//        dy = pntb1 -> y - pnta2 -> y;
//        dd = dx * dx + dy * dy;
//        if(dist < dd)dist = dd;
//    
//        dist = sqrt(dist);
//    
//        if(dist > dda + ddb)return(NO);
//        return(YES);
//      }
//    
//      ta = (dbx * dy - dby * dx) / det;
//      tb = (dax * dy - day * dx) / det;
//    
//      if(ta < 0.0 || ta > 1.0)return(NO);
//      if(tb < 0.0 || tb > 1.0)return(NO);
//    
//      return(YES);
//    }
//    
//    
//    

static int project_point_to_ellipse (point_rec *point, ellipse_rec *ellipse, point_rec *pnt, double *t_par)
{
  point_rec grad, vec, point_spec, pnt_spec, p;
  double xx, yy, dist, norm, rate1 = 0.25, rate2 = 1.5, delta_t = EPSILON_T;
  double t = 0.5, dt, dtt, ddt, d_t = 0.0, max_dt = 0.25, rate;
  int sign = 0, sign_x = 0, sign_y = 0, iter = 0, quadrant = 0;

  /* move point to the first quadrant */
  transform_from_global_to_local (point, &p, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
  xx = p.x;
  yy = p.y;

  if (xx == 0.0 && yy == 0.0)error_message (GENERAL_ERROR, "Point is in center of inclusion %d", ellipse -> id);
  if (ellipse -> max != ellipse -> min) {
    if (yy == 0.0) {
      rate = ellipse -> min / ellipse -> max;
      if (fabs (xx) < ellipse -> max * (1.0 - rate * rate))
    error_message (GENERAL_ERROR, "Point lies on major axis of inclusion %d", ellipse -> id);
    }
  }

  p.x = fabs (xx);
  p.y = fabs (yy);

  if (xx != 0.0)sign_x = (xx > 0.0) ? 1 : -1;
  if (yy != 0.0)sign_y = (yy > 0.0) ? 1 : -1;

  if (sign_x > 0) {
    if (sign_y > 0)
      quadrant = 1;
    else
      quadrant = 4;
  }
  else {
    if (sign_y > 0)
      quadrant = 2;
    else
      quadrant = 3;
  }

  if (sign_y == 0)sign = sign_x;
  if (sign_x == 0)sign = sign_y;

  if (sign != 0) {
    if (sign > 0)
      quadrant = 1;
    else
      quadrant = 3;
  }

  transform_from_local_to_global (&point_spec, &p, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);

  /* make projection (in the first quadrant) */
  while (YES) {
    get_ellipse_point (ellipse, t, &pnt_spec);
    norm = sqrt (get_ellipse_gradient (ellipse, t, &grad));
    sub_vec (vec, point_spec, pnt_spec);
    dist = dot_product (vec, grad) / norm;

    dt = dist / norm;
    if (fabs (dt) < delta_t)break;

    if ((ddt = fabs (dt)) > max_dt)dt /= ddt / max_dt;

    if (iter != 0) {
      if (dt * d_t < 0.0) {
    if (fabs (dt) > rate1 * fabs (d_t))dt = -d_t * rate1;
      }
      else {
    if (fabs (dt) > rate2 * fabs (d_t))dt = d_t * rate2;
      }
    }

    dtt = d_t = dt;

    if (dt < 0.0) {
      if (t + dt < 0.0)dtt = (0.0 - t) / 2.0;
    }
    if (dt > 0.0) {
      if (t + dt > 1.0)dtt = (1.0 - t) / 2.0;
    }

    iter++;
    t += (d_t = dtt);

    if (t < 0.0)t = 0.0;
    if (t > 1.0)t = 1.0;
  }

  /* move projection to original quadrant */
  transform_from_global_to_local (&pnt_spec, &p, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
  p.x *= sign_x;
  p.y *= sign_y;
  transform_from_local_to_global (pnt, &p, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);

  //note: t is continuous (but not monotonic)
  //         t = 0 on major axis
  //         t = 1 on minor axis
  //         to prevent confusion, quadrant is returned
  //         x+,y+   1
  //         x-,y+   2
  //         x-,y-   3
  //         x+,y-   4
  //         zeros are included into quadrants 1 and 3 

  *t_par = t;
  return(quadrant);
}

//    
//    
//    
static void
get_ellipse_point_quadrant (ellipse_rec *ellipse, int quadrant, double t, point_rec *pnt)
{
  point_rec p;

  get_ellipse_point (ellipse, t, pnt);

  /* move point to quadrant */
  if (quadrant != 1) {
    transform_from_global_to_local (pnt, &p, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
    if (quadrant == 2 || quadrant == 3)p.x *= -1.0;
    if (quadrant == 4 || quadrant == 3)p.y *= -1.0;
    transform_from_local_to_global (pnt, &p, &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
  }
}
//    
//    
//    
//    static double
//    get_ellipse_gradient_quadrant(ellipse_rec *ellipse, int quadrant, double t, point_rec *grad)
//    {
//      point_rec pnt, g;
//      double norm;
//    
//      norm = get_ellipse_gradient(ellipse, t, grad);
//    
//      /* move grad to quadrant */
//      if(quadrant != 1){
//        pnt.x = pnt.y = 0.0;
//        transform_from_global_to_local(grad, &g, &pnt, ellipse -> cos_angle, ellipse -> sin_angle);
//        if(quadrant == 2 || quadrant == 3)g.x *= -1.0;
//        if(quadrant == 4 || quadrant == 3)g.y *= -1.0;
//        transform_from_local_to_global(grad, &g, &pnt, ellipse -> cos_angle, ellipse -> sin_angle);
//      }
//    
//      return(norm);
//    }
//    
//    
//    
//    static void
//    get_ellipse_normal_quadrant(ellipse_rec *ellipse, int quadrant, double t, point_rec *normal)
//    {
//      point_rec pnt, n;
//    
//      get_ellipse_normal(ellipse, t, normal);
//    
//      /* move normal to quadrant */
//      if(quadrant != 1){
//        pnt.x = pnt.y = 0.0;
//        transform_from_global_to_local(normal, &n, &pnt, ellipse -> cos_angle, ellipse -> sin_angle);
//        if(quadrant == 2 || quadrant == 3)n.x *= -1.0;
//        if(quadrant == 4 || quadrant == 3)n.y *= -1.0;
//        transform_from_local_to_global(normal, &n, &pnt, ellipse -> cos_angle, ellipse -> sin_angle);
//      }
//    }
//    
//    
//    
static void get_ellipse_point (ellipse_rec *ellipse, double t, point_rec *pnt)
{
  double weight = 0.707106781;
  double B0, B1, B2, Bx, By;
  double sum_B;

  B0 = B_0 (t);
  B1 = B_1 (t) * weight;
  B2 = B_2 (t);

  sum_B = B0 + B1 + B2;

  Bx = B0 * ellipse -> pnt0.x + B1 * ellipse -> pnt1.x + B2 * ellipse -> pnt2.x;
  By = B0 * ellipse -> pnt0.y + B1 * ellipse -> pnt1.y + B2 * ellipse -> pnt2.y;

  pnt -> x = Bx / sum_B;
  pnt -> y = By / sum_B;
}

static double get_ellipse_gradient (ellipse_rec *ellipse, double t, point_rec *grad)
{
  double weight = 0.707106781;
  double B0, B1, B2, dB0, dB1, dB2;
  double sum_B, sum_dB, sum_B2, Bx, By, dBx, dBy;

  B0 = B_0 (t);
  B1 = B_1 (t) * weight;
  B2 = B_2 (t);

  dB0 = dB_0 (t);
  dB1 = dB_1 (t) * weight;
  dB2 = dB_2 (t);

  sum_B = B0 + B1 + B2;
  sum_dB = dB0 + dB1 + dB2;

  sum_B2 = sum_B * sum_B;

  Bx = B0 * ellipse -> pnt0.x + B1 * ellipse -> pnt1.x + B2 * ellipse -> pnt2.x;
  By = B0 * ellipse -> pnt0.y + B1 * ellipse -> pnt1.y + B2 * ellipse -> pnt2.y;

  dBx = dB0 * ellipse -> pnt0.x + dB1 * ellipse -> pnt1.x + dB2 * ellipse -> pnt2.x;
  dBy = dB0 * ellipse -> pnt0.y + dB1 * ellipse -> pnt1.y + dB2 * ellipse -> pnt2.y;

  grad -> x = (dBx * sum_B - Bx * sum_dB) / sum_B2;
  grad -> y = (dBy * sum_B - By * sum_dB) / sum_B2;

  return(size_vec (*grad));
}

//    
//    
//    
//    static void
//    get_ellipse_normal(ellipse_rec *ellipse, double t, point_rec *normal)
//    {
//      point_rec grad;
//      double norm;
//    
//      norm = get_ellipse_gradient(ellipse, t, &grad);
//      norm = sqrt(norm);
//    
//      normal -> x = grad.y;
//      normal -> y = -grad.x;
//    
//      div_vec(*normal, norm);
//    }
//    
//    
//    

static void transform_from_global_to_local (point_rec *global, point_rec *local, point_rec *origin, double cos_angle, double sin_angle)
{
  double delta_x, delta_y;

  delta_x = global -> x - origin -> x;
  delta_y = global -> y - origin -> y;

  local -> x = cos_angle * delta_x + sin_angle * delta_y;
  local -> y = cos_angle * delta_y - sin_angle * delta_x;
}

void transform_from_local_to_global (point_rec *global, point_rec *local, point_rec *origin, double cos_angle, double sin_angle)
{
  global -> x = cos_angle * local -> x - sin_angle * local -> y + origin -> x;
  global -> y = cos_angle * local -> y + sin_angle * local -> x + origin -> y;
}

//    
//    
//    
//    static char *
//    get_next_record(FILE *file, char *file_name, char *err_msg)
//    {
//      char *buf = NULL;
//    
//      if((buf = fgets(line_buffer, BUFFER_SIZE, file)) == NULL){
//        if(!feof(file))error_message(FILE_READ_ERROR, "File %s reading error", file_name);
//        if(err_msg == NULL)return(NULL);
//        error_message(INPUT_ERROR, "%s", err_msg);
//      }
//      return(buf);
//    }
//    
//    
//    
//    static char *
//    get_next_relevant_record(FILE *file, char *file_name, char *err_msg)
//    {
//      char *buf = NULL, *character = NULL;
//      
//      do{
//        if((buf = get_next_record(file, file_name, err_msg)) == NULL)return(NULL);
//        character = line_buffer;
//        while(*character == SPACE || *character == TABULATOR)character++;
//      }while(*character == NEWLINE || *character == REMARK);
//      return(buf);
//    }
//    
//    
//    
//    static void
//    write_record(FILE *file, char *file_name, char *buffer)
//    {
//      if(fprintf(file, "%s", buffer) < 0)
//        error_message(FILE_WRITE_ERROR, "File %s writing error", file_name);
//    }
//    
//    
//    


// Lukasovo - odstraneno static
void error_message (int exit_code, const char *format, ...)
{
  char buffer[1024];
  va_list args;

  va_start (args, format);
  vsprintf (buffer, format, args);
  va_end (args);

  if (exit_code != WARNING && exit_code != NO_ERROR) {
    fprintf (stderr, "\nError: %s\n\n", buffer);

#ifdef ELIXIR
    if(ElixirInitialized() == YES){
      strcat(buffer, "\n\n");
      ESIEventLoop(YES, buffer);
    }
#endif

    exit (exit_code);
  }

  if (exit_code == WARNING)
    fprintf (stderr, "Warning: %s\n", buffer);
  else
    fprintf (stderr, "%s\n", buffer);
}

//    
//    
//    
//    //#ifdef ELIXIR
//    //
//    //static void
//    //draw_point(point_rec *pnt, EPixel color)
//    //{
//    //  WCRec point[1];
//    //  GraphicObj *gr_obj = NULL;
//    //  unsigned long mask = 0L;
//    //    
//    //  mask |= LAYER_MASK;
//    //  mask |= COLOR_MASK;
//    //  mask |= MTYPE_MASK;
//    //  mask |= MSIZE_MASK;
//    //
//    //  EASValsSetLayer(0);
//    //  EASValsSetColor(color);
//    //  EASValsSetMType(FILLED_CIRCLE_MARKER);
//    //  EASValsSetMSize(6);   
//    //
//    //  point[0].z = 0.0;
//    //
//    //  copy_vec(point[0], *pnt);
//    //
//    //  if((gr_obj = CreateMarker3D(point)) == NULL){
//    //    ERptErrMessage(ELIXIR_ERROR_CLASS, 2, "Failure of graphics creation", ERROR_GRADE);
//    //    exit(ESISetExitCode(2));
//    //  }
//    //
//    //  EGWithMaskChangeAttributes(mask, gr_obj);
//    //  EMAddGraphicsToModel(ESIModel(), gr_obj);
//    //}
//    //
//    //
//    //
//    //static void
//    //draw_line(point_rec *pnt1, point_rec *pnt2, EPixel color)
//    //{
//    //  WCRec point[2];
//    //  GraphicObj *gr_obj = NULL;
//    //  unsigned long mask = 0L;
//    //    
//    //  mask |= LAYER_MASK;
//    //  mask |= COLOR_MASK;
//    //  mask |= STYLE_MASK;
//    //  mask |= SHRINK_MASK;
//    //
//    //  EASValsSetLayer(0);
//    //  EASValsSetColor(color);
//    //  EASValsSetFillStyle(SOLID_STYLE);
//    //  EASValsSetShrink(1.0);
//    //
//    //  point[0].z = point[1].z = 0.0;
//    //
//    //  copy_vec(point[0], *pnt1);
//    //  copy_vec(point[1], *pnt2);
//    //
//    //  if((gr_obj = CreateLine3D(point)) == NULL){
//    //    ERptErrMessage(ELIXIR_ERROR_CLASS, 2, "Failure of graphics creation", ERROR_GRADE);
//    //    exit(ESISetExitCode(2));
//    //  }
//    //        
//    //  EGWithMaskChangeAttributes(mask, gr_obj);
//    //  EMAddGraphicsToModel(ESIModel(), gr_obj);
//    //}
//    //
//    //
//    //
//    //
//    //static void
//    //draw_elliptic_arc(point_rec *pnt1, point_rec *pnt2, point_rec *pnt3, double w1, double w2, double w3, EPixel color)
//    //{
//    //  WCRec point[3];
//    //  float weight[3];
//    //  GraphicObj *gr_obj = NULL;
//    //  unsigned long mask = 0L;
//    //    
//    //  mask |= LAYER_MASK;
//    //  mask |= COLOR_MASK;
//    //  mask |= STYLE_MASK;
//    //  mask |= WIDTH_MASK;
//    //  mask |= SHOW_POLY_MASK;
//    //  mask |= TESSEL_INTERVALS_MASK;
//    //
//    //  EASValsSetLayer(0);
//    //  EASValsSetColor(color);
//    //  EASValsSetFillStyle(SOLID_STYLE);
//    //  EASValsSetLineWidth(0);
//    //  EASValsSetShowPoly(NO);
//    //  EASValsSetTesselIntervals(32);
//    //
//    //  point[0].z = point[1].z = point[2].z = 0.0;
//    //
//    //  copy_vec(point[0], *pnt1);
//    //  copy_vec(point[1], *pnt2);
//    //  copy_vec(point[2], *pnt3);
//    //
//    //  weight[0] = w1;
//    //  weight[1] = w2;
//    //  weight[2] = w3;
//    //
//    //  if((gr_obj = CreateRBezC3D(3, point, weight)) == NULL){
//    //    ERptErrMessage(ELIXIR_ERROR_CLASS, 2, "Failure of graphics creation", ERROR_GRADE);
//    //    exit(ESISetExitCode(2));
//    //  }
//    //        
//    //  EGWithMaskChangeAttributes(mask, gr_obj);
//    //  EMAddGraphicsToModel(ESIModel(), gr_obj);
//    //}
//    //
//    //
//    //
//    //
//    //static void
//    //draw_cylinder(cylinder_rec *cylinder, EPixel color)
//    //{
//    //  draw_line(&(cylinder -> lower_left), &(cylinder -> lower_right), color);
//    //  draw_line(&(cylinder -> upper_left), &(cylinder -> upper_right), color);
//    //  draw_line(&(cylinder -> lower_left), &(cylinder -> upper_left), color);
//    //  draw_line(&(cylinder -> lower_right), &(cylinder -> upper_right), color);
//    //}
//    //
//    //
//    //
//    //static void
//    //draw_ellipse(ellipse_rec *ellipse, EPixel color)
//    //{
//    //  point_rec point_loc[3], point_glob[3];
//    //  int i, signx[4] = {1, -1, -1, 1}, signy[4] = {1, 1, -1, -1};
//    //
//    //  point_loc[2].x = point_loc[0].y = 0.0;
//    //
//    //  for(i = 0; i < 4; i++){
//    //    point_loc[0].x = point_loc[1].x = ellipse -> max * signx[i];
//    //    point_loc[1].y = point_loc[2].y = ellipse -> min * signy[i];
//    //
//    //    transform_from_local_to_global(&(point_glob[0]), &(point_loc[0]), &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//    //    transform_from_local_to_global(&(point_glob[1]), &(point_loc[1]), &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//    //    transform_from_local_to_global(&(point_glob[2]), &(point_loc[2]), &(ellipse -> center), ellipse -> cos_angle, ellipse -> sin_angle);
//    //
//    //    draw_elliptic_arc(&(point_glob[0]), &(point_glob[1]), &(point_glob[2]), 1.0, 0.707106781, 1.0, color);
//    //  }
//    //}
//    //
//    //
//    //
//    //static void
//    //draw_crack(crack_rec *crack, EPixel color)
//    //{
//    //  draw_line(&(crack -> center), &(crack -> tip), color);
//    //}
//    //
//    //
//    //
//    //void
//    //ESICustomize(Widget parent_pane)
//    //{
//    //}
//    //
//    //#endif

} // end of namespace meso2d