/* ../xa/msh2d.c
//  MSH2D       mesh (tesselate) 2D contours functions     2017-06-20      RF
//    In:    boundaries, tolerances, 2D
//    Out:   faces (triangles), 2D
 *
 *
 * Copyright (C) 2015 CADCAM-Services Franz Reiter (franz.reiter@cadcam.co.at)
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU 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 General Public License for more details.
 *
 *
-----------------------------------------------------
TODO:
- make testmode (eg with AP_deb_stat(); using "printd" into <tempDir>/debug.dat)

-----------------------------------------------------
Modifications:
  ..
?
-----------------------------------------------------
*/
#ifdef globTag
void MSH2D(){}
#endif
/*!
\file  ../ut/msh2d.c
\brief functions for 2D-tesselation
\code

=====================================================
List_functions_start:

--------------------- USER-FUNCTIONS
MSH2D_bnd_3D_add_dbo  add DB-curve as indexed-points-curve to TessStru.pa3
MSH2D_bnd_add__       add links to ts1->ipa and boundary-record to ts1->tab
MSH2D_bnd_2D_3D       change all 3D-boundaryPoints -> 2D-boundaryPoints
MSH2D_tess_pln        define gridbox for planar surface
MSH2D_tess_grid       define gridbox from (uMin, uMax, vMin, vMax, uNr, vNr)
MSH2D_tess__          tesselate
MSH2D_pgb_get         returns list of all 2D-gridpoints
MSH2D_bp_get          get backplane for outer-boundary

--------------------- INTERNAL-FUNCTIONS
MSH2D_memspc_init     init tess - memspace
MSH2D_memspc_exit     free tess - memspace

MSH2D_init__          check if user-defined OB exists, get box for OB;
MSH2D_init_qbx        set qbx = gbx reduced (bounding OB)
MSH2D_spc_ck1         check space for additional fac's, fnb's 
MSH2D_bnd__           create all boundaries ..
MSH2D_brkl__          make breakline of boundary-segment
MSH2D_brkl_bnd__      find all edges crossing line ips-ipe
MSH2D_brkl_bnd_1_     find first edge crossing line ips-ipe
MSH2D_brkl_bnd_11     edge if1,en1 is parallel ips-ipe
MSH2D_brkl_bnd_2_     find all edges crossing line ips-ipe
MSH2D_brkl_bnd_23     faces only right of ips-ipe
MSH2D_brkl_bnd_24     faces only left of ips-ipe
MSH2D_brkl_bnd_add1   add edge to right-contour; set brkls->cp1e
MSH2D_brkl_bnd_add2   add edge to left-contour; set brkls->cp2s
MSH2D_brkl_bnd_addf   add to facelist (uniq)
MSH2D_brkl_flip       flip common edge of 2faces
MSH2D_brkl_ipba       create ipba = list of boundaryPoints
MSH2D_brkl_cpy        add new faces to mfo (overwriting deleted faces)
MSH2D_brkl_fnb_eo     get outer-neighbours-of-boundaries in brkls
MSH2D_brkl_fnb__      update fnb of new-faces and of neighbours of new-faces
MSH2D_brkl_set        set edgType of edge (fnb of the common line)
MSH2D_fdel_bnd        delete all faces inside closed contour
MSH2D_facNrMax        get max nr of faces
MSH2D_fac_del__       remove faces inside closed contour
MSH2D_fac_del_m1      mark faces deleted along closed breakline
MSH2D_fac_del_del     remove deleted faces, update fnb
MSH2D_2if_quad        get faceNrs of the 2 existing faces of quad iqx.
MSH2D_2fac_quad_b     split quad into 2 faces
MSH2D_4fac_quad_m     split quad into 4 faces
MSH2D_quad_ck1        get quadNr, column, row, if in or on gridline or gridpoint
MSH2D_quad_ck2        get the 4 quadpoints for quad itx,ity
MSH2D_quad_ckv        ck if 1|2 quads modified; set qma and faces modified; 
MSH2D_quads_p         add all gridpoints to pta and ipa
MSH2D_quads_f         mesh all quads; fix fnb
MSH2D_quad__          tesselate point in quad
MSH2D_quadNr          get nr of quads (ix * iy; 0 is first)                INLINE
MSH2D_quadInd         get quad-index (itx + (ity * MSH2D_qbx.ix)           INLINE
MSH2D_quadX           get index column (itx)                               INLINE
MSH2D_quadY           get index row (ity)                                  INLINE
MSH2D_quadLeft        get index of left-quad                               INLINE
MSH2D_quadBelow       get index of quad below                              INLINE
MSH2D_quad_is_unmod   test quad; 0=modified, 1=unmodified
MSH2D_quad_set_mf     set quadfaces modified; find edge in quad
MSH2D_quad_set_m1     set quadfaces modified
MSH2D_quad_set_mod    set quad and its faces to modified
MSH2D_quad_nf_mod     set all quads of faces to delete to modified
MSH2D_opt__           optimize - make faces less acute
MSH2D_opt_1           optimize 2 faces with common edge
MSH2D_ctess_e         tesselate contour from edges
MSH2D_ctess__         tesselate contour from polygon
MSH2D_ctess_4pt       check if contour is convex | concave (4 indexed points)
MSH2D_del__           delete unused gridpoints
MSH2D_pt__            find ipx in faces, mesh point
MSH2D_pt_gb_          find ipx on hor|vert boundary-edge, mesh point
MSH2D_pt_gbhv         check if ipx is on horizontal|vertical edge
MSH2D_pst__           alloc & init MSH2D_psta = status of the grid-points
MSH2D_pst_set         set MSH2D_psta (grid-point-status)                 INLINE
MSH2D_pst_fac_del     mark gridPoints of face i1 unused

MSH2D_view_nfac       display faces
MSH2D_view_cv_inpt    display curve of n indexed 2D-points
MSH2D_view_inpt       display n indexed 2D-points
MSH2D_dump__          UFA_fnb_dump__ & UFA_nfb_ck__
MSH2D_dump_qma        dump qma
MSH2D_dump_brkls      dump crossing-edges

MSH2D_test_ctess__    test tesselator for closed 2D-contours
MSH2D_test_qbx_disp   display gridlines
MSH2D_test_inp_CCW    test inputContours if CCW
MSH2D_test__          -

List_functions_end:
=====================================================

\endcode *//*----------------------------------------


2D - tesselation with outer-boundary, holes and a rectangular grid.

Input:
- outer-boundary and inner boundaries as indexed 2D-points, CCW
- tesselation-grid
Output:
- faces  new (grid-)points

OB = outer-boundary; CCW
IB = inner-boundary; CCW

Testmodels:
test_sur_tess1   (2 Circs ineinander)
test_msh2d_1
test_msh2d_2     (OB 24 touching support-surf)

Testprog.:
see ../xa/test_MSH2D.c
*/


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <GL/gl.h>                     // GL_TRIANGLES ..

#include "../ut/ut_geo.h"              // Point ...
#include "../ut/ut_geo_const.h"        // UT3D_CCV_NUL
#include "../ut/ut_ox_base.h"          // OGX_SET_INDEX
#include "../ut/ut_memTab.h"           // MemTab_..
#include "../ut/ut_itmsh.h"            // MSH_EDGLN_.. typedef_MemTab..
#include "../ut/ut_face.h"             // UFA_..
#include "../ut/ut_os.h"               // OS_get_bas_dir ..
#include "../ut/msh2d.h"               // TessStru

#include "../db/ut_DB.h"               // DB_get_..

#include "../xa/xa_mem.h"              // memspc51, mem_cbuf1
#include "../xa/xa_msg.h"              // MSG_*
// #include "../xa/xa_ato.h"              // ATO_getSpc_tmp__


// Externals aus ../xa/xa.c:
// extern char      AP_mod_fnam[128];
// extern Plane     WC_sur_act;            // die aktive Plane
extern int       WC_modact_ind;         // -1=primary Model is active;
                                        // else subModel is being created



/* ../ut/msh2d.h
typedef struct TessStru

//----------------------------------------------------------------
// Boundaries: tab + ipa + pa2;
//   tab        one IndTab-record in .tab = one boundary;
//   tab.ibeg   index of first point in ipa; 
//   tab.inr    nr of points of contour; last point must be same as first
//              last pointindex is (tab.ibeg + tab.inr - 1)
//   tab.typi   MSH_EDGLN_OB MSH_EDGLN_IB MSH_GRIDBOX
//   ipa        index into pta
//
//----------------------------------------------------------------
// Faces: tab + ipa + pa2 + fac
//   tab        one IndTab-record in .tab = one patch;
//   tab.ibeg   index of first point of patch in ipa;
//   tab.inr    nr of points of patch;
//              last pointindex is (tab.ibeg + tab.inr - 1)
//   tab.typi   MSH2D_PATCH  // indexed-openGl-patch
//   tab.typd   Typ_PLN (planar) or (0=Typ_Error) not-planar
//   tab.aux    GL_TRIANGLES|GL_TRIANGLE_STRIP|GL_TRIANGLE_FAN  (4-6)
//   
//----------------------------------------------------------------
// Faces (Fac3*) fac + pa2
//   i1, i2, i3 index to points (into pa2)
//   st        0-n   grid-quad-Nr (unmodified (half-)quad-face)
//             -1    face is deleted (outside OB, inside IB)
//             <MSH2D_qNr>   face is modified (not quad).
//
//----------------------------------------------------------------
// Neighbour-faces (Fac3*) fnb parallel to faces
//   i1         ind to nb-face from face-point i1 to face-point i2
//   i2         ind to nb-face from face-point i2 to face-point i3
//   i3         ind to nb-face from face-point i3 to face-point i1
//      >=0   index to neighbour-face
//       -1 = this edge is unresolved or outside or inside hole
//       -3 = this edge is a inner breakline (-MSH_EDGLN_IB)
//       -4 = this edge is a outer breakline (-MSH_EDGLN_OB)
//   st   0     all edges OK
//       16   one or more links in fnb are NOT valid
//----------------------------------------------------------------


// #ifdef _MSC_VER




#else
  // C99
static TessStru MSH2D_ts = {
                   .gbx = _GRIDBOX_NUL,
                   .tab = _MEMTAB_NUL,     // boundaries, tesselated patches
                   .ipa = _MEMTAB_NUL,     // index to points
                   .pa2 = _MEMTAB_NUL,     // 2D-points
                   .fac = _MEMTAB_NUL,     // indexed faces
                   // .pa3 = _MEMTAB_NUL,     // 3D-points
                   // .vc3 = _MEMTAB_NUL,     // normalVectors
                   // .vcz = _UT3D_VECTOR_NUL, // normalvector - only for planar
                   // .col = _ColRGB_NUL       // color|texture
                           };
#endif
*/


static TessStru     *MSH2D_ts;
static Fac3         *MFAC__ = NULL, *MFNB__ = NULL;

// status GridBoxPoints; 0=unused, 1=OK
static MemTab(char) MSH2D_mpst = {NULL, 0, 0, sizeof(char), Typ_Int1, 10, 0, 0};
//                                data rMax rNr rSiz        typ incSiz use spcTyp
static int          MSH2D_pgb_ips;        // start of gridpoints in MSH2D_ts->ipa
                                          // -1 = no gridpoints.

static double       MSH2D_tol, MSH2D_2TOL;
static double       MSH2D_eLen, MSH2D_mLen;

static GridBox      MSH2D_qbx;            // reduced GridBox (OB boxed)
static int          MSH2D_iqxMax;         // x-index of rightmost quad
static int          MSH2D_iqyMax;         // y-index of topmost quad


static int          MSH2D_qNr;            // nr of quads / fac.st: face modif.
static int          MSH2D_bnd_nr;         // nr of boundaries (OB + IB)
static int          MSH2D_it_OB;          // tab-index of OB
static char         MSH2D_is_delOB;       // 1=faces outside OB exist; 0=not
static char         *MSH2D_qma;

// static Point *MSH2D_pa;


// brklStru: table of edges crossing a breakline.

// ips,ipe       index startpoint / endpoint of breakline
// ps, pe        startpoint / endpoint of breakline
// vcl           vector ps-pe
// cf1a[c1Nr]    face-boundary on right side from ips to ipe
// ce1a[c1Nr]    edgeNr cf1a
// cf2a[c2Nr]    face-boundary on left side - from ipe to ips
// ce2a[c2Nr]    edgeNr cf2a
// fia[fiNr]     faces inside boundary
// siz           size of cf1a, ce1a, cf2a, ce2a, fia
// ixf,ixe       ixTyp=0: last edge crossing ips-ipe (cp1e-cp2s)
//               ixTyp=1: last edge right of ips-ipe (cp1e-ipxf)
// ixlf,ixle     ixTyp=1: last edge left of ips-ipe (cp2s-ipxf)
// ixTyp         type of intersection;
//                 0 = normal crossing (intersection)
//                 1 = edge parallel ips-ipe; face is right of line;
//                 2 = edge parallel ips-ipe; face is left of line;
//                 3 = ips-ipe on outer-boundary; faces right of line;
//                 4 = ips-ipe on outer-boundary; faces left of line;
//                16 = finished, OK.
// cp1e          endPt of right-contour
// cp2s          startPt of left-contour
// ipxf          ixTyp=1: point of ixf/ixlf on ips-ipe
typedef struct { int c1Nr, c2Nr, fiNr, ips, ipe, siz, ixTyp, edgTyp;
                 int *cf1a, *cf2a, *fia, ixf, ixe, ixlf, ixle, cp1e, cp2s, ipxf;
                 char *ce1a, *ce2a;
                 Point2 *ps, *pe, *pta;
                 Vector2  vcl;
               } brklStru;




// PROTOTYPES:
  double MSH2D_edgLen (int ip1, int ip2);
  int MSH2D_brkl_bnd_add2 (brklStru *brkls, int if1, int en1);
  int MSH2D_brkl_bnd_addf (brklStru *brkls, int if1);


// MSH2D_quad_is_unmod   test quad; 0=modified, 1=unmodified
#define MSH2D_quad_is_unmod(iqx) MSH2D_qma[iqx]



//================================================================
  int MSH2D_bp_get (Vector *vz, TessStru *ts1) {
//================================================================
// MSH2D_bp_get       get normalVector and backplane for outer-boundary
/// backplane: see UT3D_bp_vcz

  int      itNr, *ipa, *ia, ib, ptNr, bp, ip1, iia[6];
  Point    *pta3, *pa, pb1, pb2;
  IndTab   *itTab, *itAct;


  printf("MSH2D_bp_get \n");


  itNr  = MEMTAB_IND (&ts1->tab);  // nr of contours
  itTab = MEMTAB_DAT (&ts1->tab);  // contours
  pta3 = MEMTAB_DAT (&ts1->pa3);
  ipa  = MEMTAB_DAT (&ts1->ipa);   // indices into pa3 & pa2 & vc3


  // get points of outer-boundary
  // loop tru all bounds in .tab
  for(ib=0; ib<itNr; ++ib) {
    itAct = &itTab[ib];
    if(itAct->typi != MSH_EDGLN_OB) continue;
      // printf(" nxt-bnd ib=%d\n",ib);

    // ia = &ipa[itAct->ibeg];       // first index
    ip1 = ipa[itAct->ibeg];
    pa = &pta3[ip1];
    ptNr = itAct->iNr;            // nr of points
    goto L_1;
  }

  return -1;


  L_1:


// using pta3 not indexed !
  // get 3D-box from npt
  UT3D_box_ix_npt (&pb1, &pb2, iia, ptNr, pa);
    UT3D_stru_dump (Typ_PT, &pb1, " &pb1");
    UT3D_stru_dump (Typ_PT, &pb2, " &pb2");

  // get backPlane from 3D-box
  bp = UT3D_bp_2pt (&pb1, &pb2);
    printf(" bp=%d\n",bp);

  // get normalvector plb->vz from backPlane and extreme-points
  UT3D_vc_perp_npt_bp_ix (vz, pa, iia, bp);
    UT3D_stru_dump (Typ_VC, vz, " vz");


  // get backplane of points
  // bp =  UT3D_bp_inpt (ptNr, ia, pta3);

    printf("ex-MSH2D_bp_get bP=%d\n",bp);

  return bp;

}


//=====================================================================
   int MSH2D_bnd_2D_3D (TessStru *ts1, int (*tr_2D_3D)()) {
//=====================================================================
// MSH2D_bnd_2D_3D       change all 3D-boundaryPoints -> 2D-boundaryPoints
// use tr_2D_3D to fill ts1->pa2 from ts1->pa3

  int      irc, i1, i2, ib, itNr, ptNr, *ipa, *ia, ipx;
  long     l1;
  IndTab   *itTab, *itAct;
  Point2   *pta2;
  Point    *pta3;



  printf("MSH2D_bnd_2D_3D \n");
  // MemTab_dump (&ts1->pa3, ".pa3-bnd_2D_3D");
  // UT3D_stru_dump (Typ_MemTab, &ts1->pa2, ".pa2-bnd_2D_3D");


  //----------------------------------------------------------------
  // get space for enough points in pta2
  ptNr = MEMTAB_IND (&ts1->pa3);
  i1 = MEMTAB_RMAX(&ts1->pa2);
  if(i1 < ptNr) {
    i2 = UTI_round_32up (ptNr - i1);
      // printf(" MemTab_add-pa2 %d\n",i2);
    MemTab_add (&ts1->pa2, &l1, NULL, i2, 2);
  }

  // set ptNr
  MEMTAB_IND(&ts1->pa2) = ptNr;


  //----------------------------------------------------------------
  itNr  = MEMTAB_IND (&ts1->tab);  // nr of contours
  itTab = MEMTAB_DAT (&ts1->tab); // contours
  pta3 = MEMTAB_DAT (&ts1->pa3);
  pta2 = MEMTAB_DAT (&ts1->pa2); // must be after MemTab_add !
  ipa  = MEMTAB_DAT (&ts1->ipa);   // indices into pa3 & pa2 & vc3


  // loop tru all bounds in .tab
  for(ib=0; ib<itNr; ++ib) {
    itAct = &itTab[ib];
    if((itAct->typi != MSH_EDGLN_OB) &&
       (itAct->typi != MSH_EDGLN_IB)) continue;
      // printf(" nxt-bnd ib=%d\n",ib);

    ia = &ipa[itAct->ibeg];       // first index
    ptNr = itAct->iNr;            // nr of points

    // loop tru all points of itAct; get 2D-point from 3D-point
    // transform points
    for(i1=0; i1<ptNr; ++i1) {
      ipx = ia[i1];
        // UT3D_stru_dump (Typ_PT, &pta3[ipx], "3D-p%d",ipx);
      tr_2D_3D (&pta2[ipx], &pta3[ipx]);
        // UT3D_stru_dump (Typ_PT2, &pta2[ipx], "2D-p%d",ipx);
    }


      // TESTBLOCK
      MSH2D_view_inpt (ptNr, ia, pta2);
      // END TESTBLOCK

  }



  return 0;

}


//===========================================================================
  int MSH2D_bnd_3D_add_dbo (TessStru *ts1, int dbTyp, long dbi, int mshTyp) {
//===========================================================================
// MSH2D_bnd_3D_add_dbo   add DB-curve as indexed-points-curve to TessStru.pa3
// Input:
//   dbTyp     DB-typ (Typ_CV(S)|Typ_CI(C) ..
//   dbi       dbi of curve to add
//   mshTyp    type of boundary - MSH_EDGLN_OB|MSH_EDGLN_IB
// Output:
//   ts1
//
    
  int     irc, ptNr, ips, iis;
  long    ld;
  Point   *pa3;
    

  // get points of obj into pa3
  irc = PRCV_npt_dbo__ (&pa3, &ptNr, dbTyp, dbi, WC_modact_ind);
  if(irc < 0) {TX_Error("MSH2D_bnd_3D_add_dbo E1");return -1;}
    // printf(" _npt_dbo__ %d %d\n",irc,ptNr);
    // GR_Disp_pTab (ptNr, pa3, SYM_STAR_S, 2);


    // TESTBLOCK
    // test sr (only for z=0,0,1
    // printf(" sr=%d\n", UT3D_sr_polc (ptNr, pa3,
               // (Vector*)&UT3D_VECTOR_Z, BCKPLN_XY));
    // END TESTBLOCK


  // remove last point if (last-point == first-point)
  if(UT3D_comp2pt(&pa3[0], &pa3[ptNr-1], UT_TOL_cv))
    --ptNr; // yes, equal -remove last


  // get startIndex for the points
  ips = MEMTAB_IND (&ts1->pa3);
    // printf("MSH2D_bnd_3D_add_dbo ptNr=%d ips=%d\n",ptNr,ips);


  // add ptNr points in pa3 to ts1->pa3
  irc = MemTab_add (&ts1->pa3, &ld, pa3, ptNr, 0);
  if(irc < 0) return irc;


  // add ptNr points starting at ips as boundary-record to ts1->tab
  // and add all links to ts1->ipa
  // and add boundary-record to ts1->tab (typ=mshTyp)
  irc = MSH2D_bnd_add__ (ts1, ips, ptNr, mshTyp);


    // TESTBLOCK
    // MemTab_dump (&ts1->tab, ".tab-0");
    // MemTab_dump (&ts1->ipa, ".ipa-0");
    // MemTab_dump (&ts1->pa3, ".pa3-0");
    // END TESTBLOCK


  return irc;

}


//================================================================
  int MSH2D_bnd_add__ (TessStru *ts1, int ips, int ptNr, int mshTyp) {
//================================================================
// MSH2D_bnd_add__     add links to ts1->ipa and boundary-record to ts1->tab;
// points must have been added before.
// Input:
//   ips       index first point in pa2
//   ptNr      nr of points
//   mshTyp    type of boundary - MSH_EDGLN_OB|MSH_EDGLN_IB
// Output:
//   ts1
//   retCod    index first link in ipa
//
// see MSHIG_init_bnd_add

  int       irc, iis;
  char      c1;
  Point     *pta;


  // printf("MSH2D_bnd_add__ ptNr=%d ips=%d\n",ptNr,ips);


/*
  // remove last point if (last-point == first-point)
  pta = MEMTAB_DAT (&ts1->pa2);
  if(MSH2D_comp2pt((Point2*)&pta[ips], (Point2*)&pta[ips+ptNr-1]))
    --ptNr; // yes, equal -remove last
*/

  // add all links into ipa
  iis = MemTabI_add_s (&ts1->ipa, ptNr, ips);
  if(iis < 0) return -1;

  // add boundary-record -> ts1->tab
  irc = MemTabIT_add__ (&ts1->tab, iis, ptNr, mshTyp);
  if(irc < 0) return irc;

    // MemTab_dump (&ts1->tab, "tab-2");
    // MemTab_dump (&ts1->ipa, "ipa-2");
    // MemTab_dump (&ts1->pa2, "pa2-2");

  return iis;

}


//================================================================
  int MSH2D_tess_pln (GridBox *gbx) {
//================================================================
// MSH2D_tess_pln     define gridbox for planar surface
// see MSH2D_tess_grid

  gbx->ix = 0;
  gbx->iy = 0;
  gbx->iz = 0;

  return 0;

}


//================================================================
  int MSH2D_tess_grid (GridBox *gbx,
                       double uMin, double uMax,
                       double vMin, double vMax, int uNr, int vNr) {
//================================================================
// define gridbox from (uMin, uMax, vMin, vMax, uNr, vNr)
// uMin, uMax = the min, max 2D-x-coords
// vMin, vMax = the min, max 2D-y-coords
// uNr,  vNr  = nr of columns / rows
// all values from support-surface
// used in MSH2D_init__
// for unlimited planar support-surface use MSH2D_tess_pln


  printf("MSH2D_tess_grid uNr=%d vNr=%d\n",uNr,vNr);


  UT3D_pt_3db (&gbx->pMin, uMin, vMin, 0.);

  gbx->ix = uNr;
  gbx->iy = vNr;
  gbx->iz = 0;

  if(uNr) gbx->dx = (uMax - uMin) / uNr;
  else    gbx->dx = 0.;
  if(vNr) gbx->dy = (vMax - vMin) / vNr;
  else    gbx->dy = 0.;
  gbx->dz = 0.;

    // UT3D_stru_dump (Typ_GridBox, gbx, "MSH2D_tess_grid");

  return 0;

}

 
//================================================================
  int MSH2D_memspc_init (TessStru *ts1, int ptNr) {
//================================================================

  int  fNr;

  fNr = ptNr * 4;   // see UFA_fNr_max


  printf("MSH2D_memspc_init fNr=%d\n",fNr);


  // get memspace for tab, ipa, p2a;
  MemTab_ini (&ts1->tab, sizeof(IndTab), Typ_IndTab, ptNr / 2);
  MemTab_ini (&ts1->ipa, sizeof(int), Typ_Int4, ptNr * 2);
  MemTab_ini (&ts1->pa2, sizeof(Point2), Typ_PT2, ptNr);
  MemTab_ini (&ts1->pa3, sizeof(Point), Typ_PT, ptNr);
  MemTab_ini (&ts1->vc3, sizeof(Vec3f), Typ_VC3F, ptNr);

  MemTab_ini (&ts1->fac, sizeof(Fac3), Typ_Fac3, fNr);
  MemTab_ini (&ts1->fnb, sizeof(Fac3), Typ_Fac3, fNr);

    // TESTBLOCK
    // UT3D_stru_dump (Typ_MemTab, &ts1->tab, "tab");
    // UT3D_stru_dump (Typ_MemTab, &ts1->ipa, "ipa");
    // UT3D_stru_dump (Typ_MemTab, &ts1->pa3, "pa3");
    // UT3D_stru_dump (Typ_MemTab, &ts1->pa2, "pa2");
    // UT3D_stru_dump (Typ_MemTab, &ts1->fac, "fac");
    // UT3D_stru_dump (Typ_MemTab, &ts1->fnb, "fnb");
    // UT3D_stru_dump (Typ_GridBox, &ts1->gbx, ".gbx-rdf");
    // printf(" tol=%f\n",ts1->tol);
    // END TESTBLOCK

  return 0;

}


//================================================================
  int MSH2D_memspc_exit (TessStru *ts1) {
//================================================================
// MSH2D_memspc_exit         free tess - memspace

  printf("MSH2D_memspc_exit \n");

  // free or reset memspaces:
  MemTab_free (&ts1->tab);
  MemTab_free (&ts1->ipa);
  MemTab_free (&ts1->pa2);
  MemTab_free (&ts1->pa3);

  // DO NOT FREE tessDat.fac AS LONG OPENGL IS USING IT !!
  MemTab_free (&ts1->fac);
  MemTab_free (&ts1->fnb);


  return 0;

}

 
//================================================================
  int MSH2D_tess__ (TessStru *ts1, int iDel, int iMod) {
//================================================================
// tesselate; do NOT modify boundaries
// iMod     0 = do not optimize
//          1 = optimize = make faces less acute   
// iDel     0=keep unused gridpoints, 1=delete

  int        irc, i1;



  printf("MSH2D_tess__ ptNr=%d ======================\n",ts1->pa2.rNr);
  // printf("  .tab.rNr=%d rMax=%d rSiz=%d\n",
            // ts1->tab.rNr,ts1->tab.rMax,ts1->tab.rSiz);
  // GL_view_npt (MEMTAB_DAT(&MSH2D_ts.pa2), MEMTAB_IND(&MSH2D_ts.pa2),
               // ATT_PT_HILI, "n", -1L); // display all points
  // MemTab_dump (&ts1->tab, ".tab-t");
  // MemTab_dump (&ts1->ipa, ".ipa-t");
  // MemTab_dump (&ts1->pa2, ".pa2-t");
  // UT3D_stru_dump (Typ_MemTab, &MSH2D_mpst, "mpst-t");
  // if(MSH2D_test_inp_CCW (ts1)) return -1; // test inputContours if CCW



  //----------------------------------------------------------------
  // copy pointers MSH2D_ts, MFAC__, MFNB__, MSH2D_tol ..
  // get MSH2D_it_OB  = tab-index of OB;
  // get MSH2D_bnd_nr = nr of contours;
  // get MSH2D_qbx = box for OB-points;
  // if OB exists: get qbx = box for OB-points;
  irc = MSH2D_init__ (ts1);  // -1=no-OB; else OB-index
  if(irc < -1) return irc;


  //----------------------------------------------------------------
  // add gridbox-points and tesselate gridbox (quads); fix fnb
  if(MEMTAB_IND(&MSH2D_ts->fac) == 0) {
    // add all gridpoints to pta and ipa
    MSH2D_quads_p (&MSH2D_qbx);

    // tesselate all quads; fix fnb
    MSH2D_quads_f (&MSH2D_qbx);
    // if(irc < 0) goto L_errEx;
  }


  //----------------------------------------------------------------
  // allocate space for MSH2D_psta = status of the grid-points
  // init MSH2D_pst. Planar does not have grid-points
  if(!MSH2D_IS_PLANAR) {
    irc = MSH2D_pst__ ();
    if(irc < 0) goto L_errEx;
  }


  //----------------------------------------------------------------
  // create OB, IB's
  if(MSH2D_bnd_nr) {

    // get memspc and setup qma = bitmap for quads 
    MSH2D_qNr = MSH2D_quadNr ();
    MSH2D_qma = (char*) MEM_alloc_tmp (MSH2D_qNr);
      // printf(" MSH2D_qNr=%d\n",MSH2D_qNr);

    // set qma to 1=unmodified (0=modified).
    for(i1=0; i1<MSH2D_qNr; ++i1) MSH2D_qma[i1] = 1;

    // create all boundaries ..
    irc = MSH2D_bnd__ ();
    if(irc < 0) goto L_errEx;
  }


  //----------------------------------------------------------------
  // optimize
  if(iMod) {
    irc = MSH2D_opt__ ();
    if(irc < 0) goto L_errEx;
      // TESTBLOCK
      // MSH2D_dump__ ("S5");
      // END TESTBLOCK
  }


  //----------------------------------------------------------------
  // delete unused gridpoints
  if(iDel) {
    irc = MSH2D_del__ ();
    if(irc < 0) goto L_errEx;
      // TESTBLOCK
      // MSH2D_dump__ ("S6");
      // END TESTBLOCK
  }


    // TESTBLOCK
    // MemTab_dump (&MSH2D_ts->fnb, ".fnb-ex-tess_-3_");
    // END TESTBLOCK


  //----------------------------------------------------------------
  return 0;


  //----------------------------------------------------------------
  L_errEx:
    MSH2D_dump__ ("tess__errEx");
    MSH2D_dump_qma ();          // display status of quads
    return irc;

}


//================================================================
  int MSH2D_pgb_get (int *ptNr, int **ia, char **stat) {
//================================================================
// MSH2D_pgb_get         returns list of all 2D-gridpoints
// - Planar does not have grid-points
// DO NOT free ia, stat
// Input:
//   ptNr      nr of points
//   ia        indexes of <ptNr> points into .pa2
//   stat      point used/unused; 1 = point is used; 0: point is unused;

  int      i1, itNr, *ipa;
  IndTab   *itTab;


  printf("MSH2D_pgb_get \n");
  // UT3D_stru_dump (Typ_MemTab, &MSH2D_mpst, "mpst-t");


  // Planar does not have grid-points
  if(MSH2D_IS_PLANAR) {*ptNr = 0; return 0;}

  ipa = MEMTAB_DAT (&MSH2D_ts->ipa);

  *ptNr = MEMTAB_IND (&MSH2D_mpst);
  *ia = &ipa[MSH2D_pgb_ips];
  *stat = MEMTAB_DAT (&MSH2D_mpst);


    // TESTBLOCK
/*
// TODO: use MSH2D_view_inpt
    printf(" ex-MSH2D_pgb_get ptNr=%d\n",*ptNr);
    {   // display unused gridPoint
      int ii=0; Point2 *pa2; Point *ptax; double zVal=0.01;
      // get data-block for points
      pa2 = MEMTAB_DAT (&MSH2D_ts->pa2);
      // get tempSpc for 3D-points
      ptax = (Point*) MEM_alloc_tmp (*ptNr * sizeof(Point));
      for(i1=0; i1 < *ptNr; ++i1) {
          printf(" stat[%d] = %d ind=%d\n",i1,(*stat)[i1],(*ia)[i1]);
        if((*stat)[i1] != 0) {
          UT3D_pt_ptz (&ptax[ii], &pa2[(*ia)[i1]], zVal);
          ++ii;
        }
      }
      GL_view_npt (ptax, ii, ATT_PT_GREEN, -1L);
      // GL_view_nsymB  (ptax, ii, ATT_COL_RED, SYM_STAR_S, -1L);
    }
*/
    // END TESTBLOCK

  return 0;

}

 
//================================================================
  int MSH2D_pst__ () {
//================================================================
// MSH2D_pst__           alloc & init MSH2D_psta = status of the grid-points
// init MSH2D_psta; 1 = used as gridpoint
//                  0 = unused gridpoint
// retCod ?

  int      irc, i1, i2, igb, itNr, ptNr, *ipa, *ia;
  long     ld;
  char     *psta;
  IndTab   *itTab, *itAct;


  // printf("MSH2D_pst__ \n");


  itNr  = MEMTAB_IND (&MSH2D_ts->tab);  // nr of contours
  itTab = MEMTAB_DAT (&MSH2D_ts->tab);  // contours
  ipa  = MEMTAB_DAT (&MSH2D_ts->ipa);   // indices into pa3 & pa2 & vc3


  // find index MSH_GRIDBOX - record in .tab
  for(i1=0; i1<itNr; ++i1) {
    if(itTab[i1].typi == MSH_GRIDBOX) {
      igb = i1;
      itAct = &itTab[i1];
      ptNr = itAct->iNr;
      // offset gridBoxPoints; not yet set if support-surface comes in 
      MSH2D_pgb_ips = itAct->ibeg;
      ia = &ipa[MSH2D_pgb_ips];
      goto L_1;
    }
  }
  // no MSH_GRIDBOX ..
  MSH2D_pgb_ips = -1;
  goto L_exit;


  //----------------------------------------------------------------
  L_1:
    // printf(" MSH2D_pst__-L_1: _pgb_ips=%d ptNr=%d\n",MSH2D_pgb_ips,ptNr);

  // MSH2D_pgb_ips = index first gridBoxPoint; ia = array pointIndexes;
  // ptNr = nr of gridBoxPoints.
    // for(i1=0;i1<ptNr;++i1) printf(" ipa[%d]=%d\n",i1,ia[i1]);


  // allocate space for ptNr grid-points
  i2 = MEMTAB_RMAX(&MSH2D_mpst);
  if(i2 < ptNr) {
    i1 = UTI_round_32up (ptNr) - i2;
      // printf(" _mpst-add-1 i1=%d i2=%d\n",i1,i2);
    irc = MemTab_add (&MSH2D_mpst, &ld, NULL, i1, 2);
    MSH2D_mpst.rNr = ptNr;
    if(irc < 0) return irc;
  }

  // init MSH2D_psta to 1 - point-is-used
  psta = MEMTAB_DAT (&MSH2D_mpst);
  for(i1=0;i1<ptNr;++i1) psta[i1] = 1;

  L_exit:
    // printf("ex-MSH2D_pst__-MSH2D_pgb_ips=%d\n",MSH2D_pgb_ips);
    // MemTab_dump (&MSH2D_mpst, "MSH2D_mpst");

  return 0;

}

 
//================================================================
  int MSH2D_bnd__ () {
//================================================================
// MSH2D_bnd__      create all boundaries ..
// Input qma: 1: unmodified quad; mesh point inside quad: 2 faces -> 4 faces.
//            0: modified quad; all faces have .st = MSH2D_qNr

  int      irc, i1, i2, ib, itNr, *ipa, *ia, ptNr, ipx, ip1, ifc1, en1;
  int      ipqa[4], itx, ity, iqx, iq1x, brkTyp;
  IndTab   *itTab, *itAct;



  // printf("MSH2D_bnd__ ------ \n");

  ipa = MEMTAB_DAT (&MSH2D_ts->ipa);   // indices into GLT_pta

  itNr = MEMTAB_IND (&MSH2D_ts->tab);  // nr of contours
  itTab = MEMTAB_DAT (&MSH2D_ts->tab); // contours



  // loop tru all bounds
  for(ib=0; ib<itNr; ++ib) {
    itAct = &itTab[ib];
    if((itAct->typi != MSH_EDGLN_OB) &&
       (itAct->typi != MSH_EDGLN_IB)) continue;
      // printf(" nxt-bnd ib=%d\n",ib);

    brkTyp = itAct->typi * -1;
    ia = &ipa[itAct->ibeg];       // first index
    ptNr = itAct->iNr;


    //----------------------------------------------------------------
    // loop tru all points of boundary itAct;
    for(i1=0; i1<ptNr; ++i1) {
      ipx = ia[i1];

        // TESTBLOCK
        // printf("---------- nxt-bnd-pt  i1=%d ipx=%d\n",i1,ipx);
        // if(ipx == 26) return -99;
        // END TESTBLOCK

      // find quad with point ipx;
      // test if 1 or 2 quads unmodified;
      // mesh quad or prepare for normal meshing (MSH2D_pt__)
      irc = MSH2D_quad__ (ipx);
      if(irc > 1) goto L_tess_b; // 2|3 quad modif; point is on hor|vert bnd
      if(irc > 0) goto L_tess_p; // 1 quad modif; normal ..
      if(irc < 0) return irc;    // -1 Error
      goto L_brkl;               // 0 OK, quad modified, point meshed


      //----------------------------------------------------------------
      L_tess_b:
      // find ipx on grid-boundary, mesh point, ..
      irc = MSH2D_pt_gb_ (ipx, irc);
      if(irc < 0) return irc;
      goto L_brkl;


      //----------------------------------------------------------------
      L_tess_p:
      // find ipx in faces, mesh point, ..
      irc = MSH2D_pt__ (ipx);
      if(irc < -90) return irc;


      //----------------------------------------------------------------
      L_brkl:
      // set set OB-segment as breakline
      if(i1 > 0) irc = MSH2D_brkl__ (ia[i1-1], ipx, brkTyp);
      if(irc < 0) return irc;

    }

    // last segment -> breakline
    irc = MSH2D_brkl__ (ia[ptNr-1], ia[0], brkTyp);
    if(irc < 0) return irc;


    // delete faces outside OB | inside IB
    if(itAct->typi == MSH_EDGLN_IB) {
      irc = MSH2D_fac_del__ (ia[ptNr-1], ia[0]);
    } else {
      // delete all faces outside of OB-boundary (OB can go along SB)
      i2 = ptNr - 1;
      for(i1=0; i1<i2; ++i1) {
           // printf(" ia[%d]=%d\n",i1,ia[i1]);
        irc = MSH2D_fac_del__ (ia[i1 + 1], ia[i1]);
        // if(irc < 0) return irc;
      }
      // delete faces outside last segment
      irc = MSH2D_fac_del__ (ia[0], ia[ptNr-1]);
      MSH2D_is_delOB = 0; // faces outside OB do NOT exist
    }
    // if(irc < 0) return irc;

  }

  return 0;

}


//================================================================
  int MSH2D_fac_del__ (int ip1, int ip2) {
//================================================================
// MSH2D_fac_del__       remove faces inside closed contour
// ip1-ip2 is oriented edge of contour; delete all faces left of this edge.
// see also MSH2D_tess2_delIB__ MSH2D_tess2_delOB__

  int      irc, fNr, fc1, es1;


  // printf("MSH2D_fac_del__ ip1=%d ip2=%d\n",ip1,ip2);


  fNr = MEMTAB_IND (&MSH2D_ts->fac);
  // fTab = MEMTAB_DAT (&MSH2D_ts->fac);


  // get fc1,es1 = oriented edge ip1-ip2
  UFA_ifac_ck_edgo (&fc1, &es1, ip1, ip2, MFAC__, fNr);
    // printf(" fc1=%d\n",fc1);
  if(fc1 < 0) return -1;


  // mark all faces around fc1 deleted (fac.st = -1)
  irc =  MSH2D_fac_del_m1 (fc1);
    // TESTBLOCK
    // UFA_fnb_dump__ (MFAC__,
                    // MFNB__,
                    // MEMTAB_IND(&MSH2D_ts->fac),  " _fac_del__");
    // END TESTBLOCK
  if(irc < 0) return -1;


  // delete all faces marked for deletion; update fnb;
  // MSH2D_delFac__ (&ts1->fac, MSH2D_fnb);
  return MSH2D_fac_del_del ();

}


//================================================================
  int MSH2D_fac_del_m1 (int if1) {
//================================================================
// MSH2D_fac_del_m1      mark faces deleted along closed breakline
// if1 = starting face; delete all neighbours.
// mark for deletion with fa[if1].st = 16;
// was MSH2D_tess2_del_f
// see MSH_fac_del_oppe

  int      inf1, inf2, inf3;
  // Fac3     *fnb;


  // printf("MSH2D_fac_del_m1 %d\n",if1);

  // fa = MEMTAB_DAT (&MSH2D_ts->fac);
  // fnb = MEMTAB_DAT(&MSH2D_ts->fnb);

  // delete face if1
  MFAC__[if1].st = -1;

    
  // get 3 neighbour-faces
  inf1 = MFNB__[if1].i1; 
  inf2 = MFNB__[if1].i2; 
  inf3 = MFNB__[if1].i3;
    // printf(" inf1=%d inf2=%d inf3=%d\n",inf1,inf2,inf3);


  if(inf1 < 0) goto L_inf2;
    if(MFAC__[inf1].st == -1) goto L_inf2;   // skip if already deleted
    MSH2D_fac_del_m1 (inf1);              // recurse
    // mark gridPoints of face i1 unused
    MSH2D_pst_fac_del (&MFAC__[inf1]);


  L_inf2:
    if(inf2 < 0) goto L_inf3;
    if(MFAC__[inf2].st == -1) goto L_inf3;   // skip if already deleted
    MSH2D_fac_del_m1 (inf2);              // recurse
    // mark gridPoints of face i2 unused
    MSH2D_pst_fac_del (&MFAC__[inf2]);


  L_inf3:
    if(inf3 < 0) goto L_exit;
    if(MFAC__[inf3].st == -1) goto L_exit;   // skip if already deleted
    MSH2D_fac_del_m1 (inf3);              // recurse
    // mark gridPoints of face i3 unused
    MSH2D_pst_fac_del (&MFAC__[inf3]);


  L_exit:

  return 0;

}


//================================================================
  int MSH2D_pt_gb_ (int ipx, int idir) {
//================================================================
// MSH2D_pt_gb_     find ipx on hor|vert boundary-edge, mesh point
// ipx is on outer-boundary (type -1);
// idir     2 = horizontal; 3 = vertical

  int    irc, ifc, ie, fNr;
  Fac3   *fax, *fnx;


  // printf("MSH2D_pt_gb_ ipx=%d idir=%d\n",ipx,idir);

  fNr = MEMTAB_IND(&MSH2D_ts->fac);
  // faa = MEMTAB_DAT(&MSH2D_ts->fac);
  // fna = MEMTAB_DAT(&MSH2D_ts->fnb);


  // loop tru all faces;
  for(ifc=0; ifc<fNr; ++ifc) {
    fnx = &MFNB__[ifc];

    // skip edges != -1
    if(fnx->i1 == -1) {
      ie = 1;
      irc = MSH2D_pt_gbhv (&MFAC__[ifc], ie, ipx, idir);
      if(irc >= 0) goto L_split;
    }
 
    if(fnx->i2 == -1) {
      ie = 2;
      irc = MSH2D_pt_gbhv (&MFAC__[ifc], ie, ipx, idir);
      if(irc >= 0) goto L_split;
    }

    if(fnx->i3 == -1) {
      ie = 3;
      irc = MSH2D_pt_gbhv (&MFAC__[ifc], ie, ipx, idir);
      if(irc >= 0) goto L_split;
    }

  }

  TX_Error("MSH2D_pt_gb_ E1-%d",ipx);
  return -1;


  //----------------------------------------------------------------
  L_split:
  // edge found; split.
      // printf("-pt_gb_ ipx=%d ifc=%d ie=%d\n",ipx,ifc,ie);
  UFA_2fac_fac_pt (&MSH2D_ts->fac, MFNB__, ipx, ifc, ie);


  return 0;

}


//================================================================
  int MSH2D_pt_gbhv (Fac3 *fax, int en1, int ipx, int idir) {
//================================================================
// MSH2D_pt_gbhv    check if ipx is on horizontal|vertical edge
// - test if edge en1 of fax is horizontal|vertical  (idir)
// - test if point ipx is one this edge
// Input:
//   fax,en1     edge to test
//   ipx          point - ist it on edge ?
//   idir         2 = horizontal; 3 = vertical
// Output:
//   retcod       -1=no, 0=yes


  int       ip1, ip2;
  Point2    *p1, *p2, *pta;


  // printf("MSH2D_pt_gbhv %d %d\n",en1,idir);

  // get points of edge
  UFA_ipse_fac_esn (&ip1, &ip2, fax, en1);

  pta = MEMTAB_DAT (&MSH2D_ts->pa2);



  //----------------------------------------------------------------
  if(idir != 2) goto L_vertic;
  // idir=2; edge is horizontal
  if(fabs(pta[ip1].y - pta[ip2].y) > MSH2D_2TOL) return -1;
  if(fabs(pta[ip1].y - pta[ipx].y) > MSH2D_2TOL) return -1;

  // test if edge goes right or left
  if((pta[ip2].x > pta[ip1].x)) {
    // edge goes right
    if(pta[ipx].x < pta[ip1].x) return -1;
    if(pta[ipx].x > pta[ip2].x) return -1;

  } else {
    // edge goes left
    if(pta[ipx].x > pta[ip1].x) return -1;
    if(pta[ipx].x < pta[ip2].x) return -1;
  }

  goto L_on;


  //----------------------------------------------------------------
  L_vertic:
  // idir=3; edge is vertical
  if(fabs(pta[ip1].x - pta[ip2].x) > MSH2D_2TOL) return -1;
  if(fabs(pta[ip1].x - pta[ipx].x) > MSH2D_2TOL) return -1;

  // test if edge goes up or down
  if((pta[ip2].y > pta[ip1].y)) {
    // edge goes up
    if(pta[ipx].y < pta[ip1].y) return -1;
    if(pta[ipx].y > pta[ip2].y) return -1;

  } else {
    // edge goes down
    if(pta[ipx].y > pta[ip1].y) return -1;
    if(pta[ipx].y < pta[ip2].y) return -1;
  }



  //----------------------------------------------------------------
  L_on:
  //   test if ipx is on edge
    // printf("ex-MSH2D_pt_gbhv-irc=0 ipx=%d ip1=%d ip2=%d\n",ipx,ip1,ip2);

  return 0;

}


//================================================================
  int MSH2D_pt__ (int ipx) {
//================================================================
// find ipx in faces, mesh point, ..
// was MSH2D_tess2_OB_pt_m


  int      irc, i1, fNr, if1, ie1, if2;
  Fac3     *fAct;
  Point2   *pta;


  // printf("MSH2D_pt__ %d\n",ipx);
  // MSH2D_dump__ ("MSH2D_pt__");


  // fa = MEMTAB_DAT (&MSH2D_ts->fac);
  fNr = MEMTAB_IND (&MSH2D_ts->fac);

  pta = MEMTAB_DAT (&MSH2D_ts->pa2);
  // fnb = MEMTAB_DAT (&MSH2D_ts->fnb);


  //----------------------------------------------------------------
  // find face of point pta[ipx]

  // loop tru all faces
  for(if1=0; if1<fNr; ++if1) {
     fAct = &MFAC__[if1];
     // skip all but unmodified faces
     if(fAct->st != MSH2D_qNr) continue;
       // printf(" pt_1-test if1=%d\n",if1);

     irc = UT2D_ck_pt_in_tria_tol (&pta[fAct->i1],&pta[fAct->i2],&pta[fAct->i3],
                                    &pta[ipx], &MSH2D_tol);
       // printf(" UT2D_ck_pt_in_tria_tol-irc=%d\n",irc);
     if(irc > 0) continue;     // 1        outside ..
     if(irc < -3) goto L_done; // -4|-5|-6 point alread exists; ignore point
     if(irc < 0)  goto L_mod1; // -1|-2|-3 point is on Edge of face
     goto L_mod2;               // 0        point inside face
  }


  //----------------------------------------------------------------
  // point is not inside face.
  // point <ipx> of a IB is outside the outer-boundary. NOT ALLOWED.
    TX_Error("******* MSH2D_pt__ E1-%d",ipx);
    return -99;


  //----------------------------------------------------------------
  // -1|-2|-3: point pta[ipx] is on edge ie1 of face if1. Create 2 * 2 new faces.
  L_mod1:
    ie1 = -irc;   // edge;

    // get if2 = neigbourFace
    if2 = UFA_if_getNf_ife (if1, ie1, MFNB__);
      // printf("L_mod1 ipx=%d if1=%d ie1=%d if2=%d fNr=%d\n",ipx,if1,ie1,if2,fNr);

    //----------------------------------------------------------------
    // if normal neighbour:
    if(if2 >= 0) {
      irc = UFA_4fac_2fac_pt (&MSH2D_ts->fac, MFNB__, ipx, if1, ie1);
      if(irc >= 0) {
        // set st of both faces to modified
        MFAC__ = MEMTAB_DAT (&MSH2D_ts->fac);
        MFAC__[if1].st = MSH2D_qNr;
        MFAC__[irc].st = MSH2D_qNr;
        MFAC__[fNr].st = MSH2D_qNr;
        MFAC__[fNr+1].st = MSH2D_qNr;
      } else {
        TX_Error("******* MSH2D_pt__ E2-%d",ipx);
        return -90;
      }

        // TESTBLOCK
        // printf(" L_mod1: irc=%d\n",irc);
        // END TESTBLOCK
      goto L_done;

    //----------------------------------------------------------------
    } else if(if2 == -1) {
      // if2 ist on boundary:

      irc = UFA_2fac_fac_pt (&MSH2D_ts->fac, MFNB__, ipx, if1, ie1);
      // set st of faces to modified
      MFAC__[if1].st = MSH2D_qNr;
      MFAC__[fNr].st = MSH2D_qNr;  // new face

      goto L_done;
    }

    TX_Error("******* MSH2D_pt__ E3-%d",ipx);
    return -90;


  //----------------------------------------------------------------
  // 0: Faces; point pta[ipx] lies inside face if1
  L_mod2:
      // printf("L_mod2 ipx=%d if1=%d fNr=%d\n",ipx,if1,fNr);
    // Create 3 new faces.
    // Face durch 3 neue Faces (die durch den Punkt gehen) ersetzen
    irc = UFA_3fac_fac_pt (&MSH2D_ts->fac, MFNB__, ipx, if1);
    // set st of both faces to modified
    if(irc >= 0) {
      MFAC__ = MEMTAB_DAT (&MSH2D_ts->fac);
      MFAC__[if1].st = MSH2D_qNr;
      MFAC__[fNr].st = MSH2D_qNr;
      MFAC__[fNr+1].st = MSH2D_qNr;
    }

      // TESTBLOCK
      // END TESTBLOCK
    goto L_done;


  //----------------------------------------------------------------
  L_done:

    // TESTBLOCK
    // printf("ex-MSH2D_pt__ irc=%d ipx=%d fNr=%d\n",
           // irc,ipx,MEMTAB_IND (&MSH2D_ts->fac));
    // MSH2D_dump__ ("ex-pt__");
    // if(ipx == 0) return -99;
    // END TESTBLOCK



  return irc;

}

//================================================================
  int MSH2D_quads_p (GridBox *qbx) {
//================================================================
// add all gridpoints to pta and ipa
// was MSH2D_tess_grd_p


  int    irc, i1, i2, ips, iis, uNr, vNr, ptNr;
  long   l1;
  double uMin, vMin, du, dv, dx, dy;
  Point2 pt1, *pta;
  IndTab it1;


  // UT3D_stru_dump (Typ_GridBox, qbx, "MSH2D_quads_p-qbx");

  uMin = qbx->pMin.x;
  vMin = qbx->pMin.y;

  uNr = qbx->ix;
  vNr = qbx->iy;

  du = qbx->dx;
  dv = qbx->dy;
    // printf("MSH2D_quads_p uMin=%f vMin=%f uNr=%d vNr=%d du=%f dv=%f\n",
                          // uMin,   vMin,   uNr,   vNr,   du,   dv);


  //================================================================
  // create all gridpoints into tempspace pta
  // see also UT3D_ptgrid_sbsp
  // create uNr*vNr gridpoints; add to pointTable;

  // get space for uNr*vNr points
  ptNr = (uNr + 1) * (vNr + 1);
  pta = (Point2*) MEM_alloc_tmp (ptNr * (int)sizeof(Point2));
  if(!pta) TX_Error("MSH2D_quads_p E1");


  iis = 0;
  dy = vMin;
  for(i1=0; i1 <= vNr; ++i1) {
    // next v = row = hor line
    dx = uMin;
    for(i2=0; i2 <= uNr; ++i2) {
      // next u = col = next vert line
      pta[iis].y = dy;
      pta[iis].x = dx;
      // pta[iis].z = 0.;
        // APT_disp_SymB2 (SYM_TRI_B, ATT_COL_BLACK, &pta[iis]);
        // UT3D_stru_dump (Typ_PT2, &pta[iis], "p-%d",iis);
      dx += du;
      ++iis;
     }
     dy += dv;
   }


  //----------------------------------------------------------------
  // get startpos in pa2
  ips = MEMTAB_IND (&MSH2D_ts->pa2);
    // printf(" ips=%d\n",ips);


  // add points
  irc = MemTab_add (&MSH2D_ts->pa2, &l1, pta, ptNr, 0);
  if(irc < 0) {TX_Error("MSH2D_quads_p EOM-pa2"); return -1;}


  // add all links to gridpoints to ipa
  // add boundary-record 'MSH_GRIDBOX' to tab
  irc = MSH2D_bnd_add__ (MSH2D_ts, ips, ptNr, MSH_GRIDBOX);
  if(irc < 0) {TX_Error("MSH2D_quads_p EOM-tab"); return -1;}
  MSH2D_pgb_ips = irc;


  //----------------------------------------------------------------
  // update pa3.rNr = pa2.rNr
  MSH2D_ts->pa3.rNr = MSH2D_ts->pa2.rNr;

  //----------------------------------------------------------------
    // MemTab_dump (&MSH2D_ts->tab, ".tab-ex-quads_p");

  return 0;

}


//================================================================
  int MSH2D_quads_f (GridBox *qbx) {
//================================================================
// mesh all quads; fix fnb
// see MSH2D_tess_grd_f MSH2D_tess_grd4_f

  int      i1, i2, ii, ix, iy, p1Nr, p2Nr, ip1, ip2, ip3, ip4, igx, ifc, itNr;
  int      ixe, iye;
  Fac3     *fnbAct;
  IndTab   *itTab;


  // printf("MSH2D_quads_f \n");

  ix = qbx->ix;
  iy = qbx->iy;

  p1Nr = iy;
  p2Nr = ix;

  igx = 0;

  ifc = MEMTAB_IND (&MSH2D_ts->fac);  // face-index
  itNr = MEMTAB_IND (&MSH2D_ts->tab);  // nr of contours
  itTab = MEMTAB_DAT (&MSH2D_ts->tab); // contours
  // fnb__ = MEMTAB_DAT (&MSH2D_ts->fnb);


  //----------------------------------------------------------------
  // get ii = startIndex in ipa
  ii = -1;
  for(i1=0; i1<itNr; ++i1) {
    if(itTab[i1].typi == MSH_GRIDBOX) ii = itTab[i1].ibeg;
  }
  if(ii < 0) {printf("**** MSH2D_quads_f E1\n"); return -2;}
    // printf(" ii=%d\n",ii);



  //----------------------------------------------------------------
  // create faces
  // i1 = loop vertical, i2 = loop horizontal
  iye = p1Nr -1;
  ixe = p2Nr -1;
  for(i1=0; i1<p1Nr; ++i1) {
    // loop horizontal with i2
    for(i2=0; i2<p2Nr; ++i2) {
      // CCW:
      ip1 = ii;
      ip2 = ii + 1;
      ip3 = ii + ix + 2;
      ip4 = ii + ix + 1;
        // printf(" iga = %d %d %d %d\n",ip1,ip2,ip3,ip4);
        // UT3D_stru_dump (Typ_PT, &pa[ip1], "ip1");
        // UT3D_stru_dump (Typ_PT, &pa[ip2], "ip2");
        // UT3D_stru_dump (Typ_PT, &pa[ip3], "ip3");
        // UT3D_stru_dump (Typ_PT, &pa[ip4], "ip4");

      // set gridboxNumber into (Fac3).st
      // igx = UTI_I32_2I16 (i1, i2);

      // add new Fac3 ip1,ip2,ip3
      UFA_add_fac_st (&MSH2D_ts->fac, ip1, ip2, ip3, igx);

      // add new Fac3 ip1,ip3,ip4
      UFA_add_fac_st (&MSH2D_ts->fac, ip1, ip3, ip4, igx);

      // create fnb for lower-right
      fnbAct = &MFNB__[ifc]; // ip1, ip2, ip3

      if(i1 == 0) fnbAct->i1 = -1;            // bottom border
      else fnbAct->i1 = ifc - (ix * 2) + 1;

      if(i2 == ixe) fnbAct->i2 = -1;          // right border
      else fnbAct->i2 = ifc + 3;

      fnbAct->i3 = ifc + 1;

      fnbAct->st = 0;

      ++ifc;
      // create fnb for upper-left
      fnbAct = &MFNB__[ifc]; // ip1, ip3, ip4

      fnbAct->i1 = ifc - 1;

      if(i1 == iye) fnbAct->i2 = -1;          // top border
      else fnbAct->i2 = ifc + (ix * 2) - 1;

      if(i2 == 0) fnbAct->i3 = -1;            // left border
      else fnbAct->i3 = ifc - 3;

      fnbAct->st = 0;                

      ++ifc;
      ++ii;
      ++igx;
    }
    ++ii;
  }


  //----------------------------------------------------------------
  if(MSH2D_qNr < 2) {
    // for 1x1 do not use qma ..
    MFAC__ = MEMTAB_DAT (&MSH2D_ts->fac);
    MFAC__[0].st = MSH2D_qNr;
    MFAC__[1].st = MSH2D_qNr;
  }


  //----------------------------------------------------------------
  // set nfb.rNr
  if(MSH2D_ts->fnb.rMax < MSH2D_ts->fac.rNr) {
    TX_Error("MSH2D_quads_f E-nfb-1");
    return -99;
  }
  MSH2D_ts->fnb.rNr = MSH2D_ts->fac.rNr;
    // printf(" MSH2D_ts->fnb.rNr=%d\n",MSH2D_ts->fnb.rNr);


  return 0;

}


//================================================================
  int MSH2D_quad__ (int ipx) {
//================================================================
// MSH2D_quad__     tesselate point in quad
// if ipx is in quad already modified return -1 (do normal mesh)
// test if ipx is inside quad or on quad-boundary or on quad-point
// Output:
//   retCod   1   use MSH2D_pt__ (quad(s) already modified or prepared for
//            2   ipx is on outer horizontal bound; quad already modified;
//            3   ipx is on outer vertical bound; quad already modified;
//            0   point ipx meshed;
//           -1   internal error ..
//                

  int    irc, itx, ity, iqx, iqm, iq1x, ipqa[4], ip1;
  int    ifc1, ifc2, en1;
  // Point2 *pta, *ptx;
  // Fac3   *fnb__;


  // pta   = MEMTAB_DAT (&MSH2D_ts->pa2);   // points
  // ptx   = &pta[ipx];
  // fnb__ = MEMTAB_DAT (&MSH2D_ts->fnb);


  // printf("MSH2D_quad__ %d\n",ipx);


  // test if ipx is in quad or on quad-boundary or on quad-point
  irc = MSH2D_quad_ck1 (&itx, &ity, ipx);
  // irc = 0=inside-quad, 1=on-vertLn, 2=on-horLn, 3=on-pt
  // ipqa = the 4 quadpoints for quad itx,ity
  if(irc < 0) return irc;


  // get iqx = quad-nr. low-left=0 ..
  iqx = MSH2D_quadInd (itx, ity);
    // printf(" _quadInd iqx=%d\n",iqx);


  //----------------------------------------------------------------
  if(irc == 0) {
    // 0 - pt-is-inside-quad: check if quad already meshed
    // ck quad-modified; 1=no, 0=yes
    if(MSH2D_qma[iqx] == 0) return 1; // quad already meshed; 
    // get ipqa = quad-points
    MSH2D_quad_ck2 (ipqa, itx, ity);
    // 1 = quad not yet meshed: mesh quad into 4 faces
    // set f.st=modif.
    irc = MSH2D_4fac_quad_m (ipqa, iqx, ipx);
    if(irc < 0) return irc;
    return 0;



  //----------------------------------------------------------------
  } else if(irc == 1) {
    // 1 = point is on vertical edge of quad itx;

    // if ipx is on leftmost-line (vertical, itx=0)
    if(itx == 0) {
      // ck quad-modified; 1=no, 0=yes
      if(!MSH2D_qma[iqx]) return 3; // quad already meshed ..
      // split face into 2 faces; ipx is on outer boundary of the grid.
      // set f.st=modif.
      irc = MSH2D_2fac_quad_b (ipqa, iqx, ipx, 3);
      if(irc < 0) return irc;
      return 0;     // point meshed ..
    }

    // if ipx is on rightmost-line (vertical)
    if(itx == MSH2D_iqxMax) {
      iq1x = iqx - 1;
      // test left-quad
      if(MSH2D_qma[iq1x] == 0) return 3; // quad is already meshed; ..
      // split face into 2 faces; ipx is on outer boundary of the grid.
      // set f.st=modif.
      irc = MSH2D_2fac_quad_b (ipqa, iq1x, ipx, 1);
      if(irc < 0) return irc;
      return 0;     // point meshed ..
    }

    // get iq1x = quadNr of quad left
    // -1 = iqx is on left boundary; invalid.
    iq1x = MSH2D_quadLeft (iqx, itx);
      // printf(" iq1x-Left = %d\n",iq1x);

    // check quads
    if((!MSH2D_qma[iqx]) && (!MSH2D_qma[iq1x])) return 1; // both quads modif.

    // one quad is unmodif; get ipqa = quad-points
    MSH2D_quad_ck2 (ipqa, itx, ity);

    // set quads and its faces to modified (qma=0, f.st=modif.)
    // and get common-edge as ifc1,en1
    irc = MSH2D_quad_ckv (&ifc1, &en1, iq1x, iqx,
                               ipqa[0]+MSH2D_pgb_ips, ipqa[3]+MSH2D_pgb_ips);

    // split edge ifc1,en1 at point ipx
    irc = UFA_4fac_2fac_pt (&MSH2D_ts->fac, MFNB__, ipx, ifc1, en1);
    if(irc < 0) return irc;
    return 0;


  //----------------------------------------------------------------
  } else if(irc == 2) {
    // 2 = point is on horizontal edge of quad;

    // if ipx is on bottom-line (horizontal, ity=0)
    if(ity == 0) {
      // ck quad-modified; 1=no, 0=yes
      if(MSH2D_qma[iqx] == 0) return 2; // quad already meshed ..
      // split face into 2 faces; ipx is on outer boundary of the grid.
      // set f.st=modif.
      irc = MSH2D_2fac_quad_b (ipqa, iqx, ipx, 0);
      if(irc < 0) return irc;
      return 0;     // point meshed ..
    }

    // if ipx is on topmost-line (horizontal)
    if(ity == MSH2D_iqyMax) {
      // test quad-below if modif.; 1=no, 0=yes
      iqm = iqx-MSH2D_iqxMax;
      if(!MSH2D_qma[iqm]) return 2; // quad is already meshed; ..
      // split face into 2 faces; ipx is on outer boundary of the grid.
      irc = MSH2D_2fac_quad_b (ipqa, iqm, ipx, 2);
      if(irc < 0) return irc;
      return 0;     // point meshed ..
    }

    // get iq1x = quadNr of quad below
    // -1 = iqx is on bottom-boundary; invalid.
    iq1x = MSH2D_quadBelow (iqx);
      // printf(" iq1x-Below = %d\n",iq1x);

    // check quads
    if((!MSH2D_qma[iqx]) && (!MSH2D_qma[iq1x])) return 1; // both quads modif.

    // one quad is unmodif; get ipqa = quad-points
    MSH2D_quad_ck2 (ipqa, itx, ity);

    // set quads and its faces to modified (qma=0, f.st=modif.)
    // and get common-edge as ifc1,en1
    irc = MSH2D_quad_ckv (&ifc1, &en1, iq1x, iqx,
                               ipqa[1]+MSH2D_pgb_ips, ipqa[0]+MSH2D_pgb_ips);

    // split edge ifc1,en1 at point ipx
    irc = UFA_4fac_2fac_pt (&MSH2D_ts->fac, MFNB__, ipx, ifc1, en1);
    if(irc < 0) return irc;
    return 0;


  //----------------------------------------------------------------
  } else if(irc == 3) {
    // 3 = point ipx == (MSH2D_pgb_ips + itx + ((MSH2D_qbx.ix + 1) * ity))
    ip1 =  MSH2D_pgb_ips + itx + ((MSH2D_qbx.ix + 1) * ity);
      // printf(" ipx=%d ip1=%d\n",ipx,ip1);
    // replace quadPoint ip1 by boundaryPoint ipx
    // mark point ip1 unused
    MSH2D_pst_set (ip1, 0);
    // keep quads unmodif.  Replace ip1 by ipx
    UFA_ipn_chg_ip (ip1, ipx,
                    MFAC__,
                    MEMTAB_IND (&MSH2D_ts->fac));
    return 0;


  //----------------------------------------------------------------
  } else {
    TX_Error("MSH2D_quad__-E1"); return -99;
  }

  irc = 0;



  //----------------------------------------------------------------
  L_exit:

    // printf("ex-MSH2D_quad__ ipx=%d irc=%d\n",ipx,irc);

  return irc;

}


//================================================================
  int MSH2D_2fac_quad_b (int *ipqa, int iqx, int ipx, int iPos) {
//================================================================
// point ipx is on boundary of quad iqx; split quad into 2 faces
// Input:
//   ipqa      the 4 quadpoints
//   ipx       point at boundary; split here
//   iPos      boundary-type; 0=bottom-edge, 1=right-vertical,
//             2=top-edge, 3=left-vertical-edge of quad
// set f.st=modif.
// set MSH2D_qma[iqx] modif.

static int ena[]={1,2,2,3};

  int     irc, ifc1, ifc2, en1, *ipa, *ipga;
  // Fac3    *fnb;


  // printf("MSH2D_2fac_quad_b iqx=%d ipx=%d iPos=%d\n",iqx,ipx,iPos);

  ipa = MEMTAB_DAT (&MSH2D_ts->ipa);
  ipga = &ipa[MSH2D_pgb_ips];

  // find faceNr 
  irc = MSH2D_2if_quad (&ifc1, &ifc2, iqx, ipga[ipqa[1]]);
  if(irc < 0) return irc;
    // printf(" ifc1=%d ifc2=%d\n",ifc1,ifc2);

  // set status of both faces to modified
  MFAC__[ifc1].st = MSH2D_qNr;
  MFAC__[ifc2].st = MSH2D_qNr;

  // ipos=0,1; use ifc1; ipos=2,3; use ifc2;
  if(iPos > 1) ifc1 = ifc2;

  // ipos=0; ifc1, en1=1;
  // ipos=1; ifc1, en1=2; 
  // ipos=2; ifc2, en1=2;
  // ipos=3; ifc2, en1=3; 
  en1 = ena[iPos];
    // printf(" 2fac_quad_-ifc1=%d en1=%d\n",ifc1,en1);

  
  // split face ifc1 at ipx
  UFA_2fac_fac_pt (&MSH2D_ts->fac, MFNB__, ipx, ifc1, en1);

  // set quad = modified
  MSH2D_qma[iqx] = 0;

    // TESTBLOCK
    // if(ipx == 0) {
      // printf("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n");
      // return -99;
    // }
    // END TESTBLOCK


  return 0;

}


//================================================================
  int MSH2D_4fac_quad_m (int *ipqa, int iqx, int ipx) {
//================================================================
// point ipx is inside quad ipa; split quad into 4 faces
// make faces 0,1,X; 1,2,X; 2,3,X; 3,0,X.
// set MSH2D_qma[iqx] modified
// Input:
//   ipqa     points (without offset)
//   iqx      quadNr
//   ipx      new point


  int   irc, i1, fNr, if1, ie1, *ipa, *ipga, ifc1, ifc2;
  int   in0, in1, in2, in3;
  Fac3  *fnbAct;


  ipa = MEMTAB_DAT (&MSH2D_ts->ipa);
  ipga = &ipa[MSH2D_pgb_ips];


  // printf("MSH2D_4fac_quad_m ipx=%d iqx=%d\n",ipx,iqx);
  // for(i1=0;i1<4;++i1) printf(" %d iq=%d ip=%d\n",i1,ipqa[i1],ipga[ipqa[i1]]);


  // fa = MEMTAB_DAT (&MSH2D_ts->fac);
  fNr = MEMTAB_IND (&MSH2D_ts->fac);  // ind of first face

  // fnb__ = MEMTAB_DAT(&MSH2D_ts->fnb);
  fnbAct = &MFNB__[fNr];           // parallel fac


  //----------------------------------------------------------------
  // get faceNrs of the 2 existing faces of quad iqx.
  irc = MSH2D_2if_quad (&ifc1, &ifc2, iqx, ipga[ipqa[1]]);
  if(irc < 0) return irc;
    // printf(" ifc1=%d ifc2=%d fNr=%d\n",ifc1,ifc2,fNr);

  // get neighbours of quad ifc1,ifc2
  // in0 = MFNB__[ifc1].i1;
  in1 = MFNB__[ifc1].i2;
  in2 = MFNB__[ifc2].i2;
  in3 = MFNB__[ifc2].i3;
    // printf(" in1=%d in2=%d in3=%d\n",in1,in2,in3);



  //----------------------------------------------------------------
  // change face ifc1 into 0,1,X;
  MFAC__[ifc1].i1 = ipga[ipqa[0]];
  MFAC__[ifc1].i2 = ipga[ipqa[1]];
  MFAC__[ifc1].i3 = ipx;
  MFAC__[ifc1].st = MSH2D_qNr;
  // set neighbours
  // MFNB__[ifc1].i1 = in0;
  MFNB__[ifc1].i2 = ifc2;
  MFNB__[ifc1].i3 = fNr + 1;


  //----------------------------------------------------------------
  // change face ifc2 into 1,2,X;
  MFAC__[ifc2].i1 = ipga[ipqa[1]];
  MFAC__[ifc2].i2 = ipga[ipqa[2]];
  MFAC__[ifc2].i3 = ipx;
  MFAC__[ifc2].st = MSH2D_qNr;
  // set neighbours
  MFNB__[ifc2].i1 = in1;
  MFNB__[ifc2].i2 = fNr;
  MFNB__[ifc2].i3 = ifc1;
  // find fnb with oriented edge 2-1
  if(in1 >= 0) {
    // set MFNB__[in1] - edge with startNr ipga[ipqa[2]] to ifc2
    UFA_fnb_set_esn (in1, ipga[ipqa[2]], ifc2, MFAC__, MFNB__);
  }


  //----------------------------------------------------------------
  // add face 2,3,X;
  // UFA_add_fac_st (&MSH2D_ts->fac, ipga[ipqa[2]], ipga[ipqa[3]], ipx, iqx);
  UFA_add_fac_st (&MSH2D_ts->fac, ipga[ipqa[2]], ipga[ipqa[3]], ipx, MSH2D_qNr);
  // set neighbours
  MFNB__[fNr].i1 = in2;
  MFNB__[fNr].i2 = fNr + 1;
  MFNB__[fNr].i3 = ifc2;
  // find fnb with oriented edge 3-2
  if(in2 >= 0) {
    // set MFNB__[in2] - edge with startNr ipga[ipqa[2]] to fNr
    UFA_fnb_set_esn (in2, ipga[ipqa[3]], fNr, MFAC__, MFNB__);
  }
  ++fNr;


  //----------------------------------------------------------------
  // add face 3,0,X;
  // UFA_add_fac_st (&MSH2D_ts->fac, ipga[ipqa[3]], ipga[ipqa[0]], ipx, iqx);
  UFA_add_fac_st (&MSH2D_ts->fac, ipga[ipqa[3]], ipga[ipqa[0]], ipx, MSH2D_qNr);
  // set neighbours
  MFNB__[fNr].i1 = in3;
  MFNB__[fNr].i2 = ifc1;
  MFNB__[fNr].i3 = fNr - 1;
  // find fnb with oriented edge 0-3
  if(in3 >= 0) {
    // set MFNB__[in3] - edge with startNr ipga[ipqa[0]] to fNr
    UFA_fnb_set_esn (in3, ipga[ipqa[0]], fNr, MFAC__, MFNB__);
  }


  //----------------------------------------------------------------
  MSH2D_qma[iqx] = 0;  // modified


  return 0;

}


//================================================================
  int MSH2D_2if_quad (int *ifc1, int *ifc2, int iqx, int ip2) {
//================================================================
// MSH2D_2if_quad        get faceNrs of the 2 existing faces of quad iqx.
// get faceNrs of the 2 existing faces of unmodified quad iqx.
// Input:
//   iqx      quadNr 
//   ip2      2. point of quad (lower right)
// Output:
//   ifc1     lower right face
//   ifc2     upper left face


  int      i1, if1, if2, fNr;
  // Fac3     *fa;


  // printf("MSH2D_2if_quad iqx=%d ip2=%d\n",iqx,ip2);

  if(MSH2D_is_delOB) {
    // faces outsideOB already exist
    *ifc1 = iqx * 2;
    *ifc2 = *ifc1 + 1;
    goto L_exit;
  }


  //----------------------------------------------------------------
  // find faces after delete faces outside OB
  // fa = MEMTAB_DAT (&MSH2D_ts->fac);
  fNr = MEMTAB_IND (&MSH2D_ts->fac);
  if1 = -1;
  if2 = -1;

  for(i1=0; i1<fNr; ++i1) {
    if(MFAC__[i1].st != iqx) continue;
    if(if1 < 0) {
      if1 = i1;
      continue;
    }
    if2 = i1;
    goto L_found;
  }

  // not found
  TX_Error("MSH2D_2if_quad E1 %d %d",iqx,ip2);
  return -1;


  L_found:
  if(MFAC__[if1].i2 == ip2) {
    *ifc1 = if1;
    *ifc2 = if2;
  } else {
    *ifc1 = if2;
    *ifc2 = if1;
  }


  L_exit:
    // printf(" ex-2if_quad ifc1=%d ifc2=%d\n",*ifc1,*ifc2);


  return 0;

}


//=================================================================
  int MSH2D_quad_ck1 (int *itx, int *ity, int ipx) {
//=================================================================
// MSH2D_quad_ck1  get quadNr, column, row, if in or on gridline or gridpoint
// Output:
//   itx      nr grid-column
//   ity      nr grid-row
//   retCod   0   point is inside quad itx,ity
//            1   point is on vertical line itx
//            2   point is on horizontal line ity
//            3   point is on gridpoint itx,ity


  int      irc, i1, ix, iy, uNr, vNr;
  long     l1;
  double   dx, dy;
  Point2   *pta, *ptx;


  pta = MEMTAB_DAT (&MSH2D_ts->pa2);   // points
  ptx = &pta[ipx];


  // printf("MSH2D_quad_ck1 %d %f %f\n",ipx,ptx->x,ptx->y);
  // if(ipx==7){
    // printf(" MSH2D_2TOL=%f\n",MSH2D_2TOL);
    // UT3D_stru_dump (Typ_GridBox, &MSH2D_qbx, "  qbx");
  // }


  // get itx=col-index, ity=row-index, dx=dist-vert, dy=dist-hor
  MSH2D_quadX (itx, &dx, ptx);
  MSH2D_quadY (ity, &dy, ptx);
    // printf(" dx=%f dy=%f\n",dx,dy);


  irc = 0;
  if(fabs(dx) < MSH2D_tol) irc += 1;
  if(fabs(dy) < MSH2D_tol) irc += 2; // on horizontal line | on gridpoint


    // printf("ex-_quad_ck1 irc=%d itx=%d ity=%d\n",irc,*itx,*ity);
    // if(ipx == 1) return -99;

  return irc;

}


//=================================================================
  int MSH2D_quad_ck2 (int *ipqa, int itx, int ity) {
//=================================================================
// ipqa = get the 4 quadpoints for quad itx,ity
// Output:
//   ipqa     index quad-points (without offset)

  int      uNr, vNr;

  uNr = MSH2D_qbx.ix + 1;
  vNr = MSH2D_qbx.iy + 1;

  ipqa[0] = uNr * ity + itx;
  ipqa[1] = ipqa[0] + 1;
  ipqa[2] = ipqa[1] + uNr;
  ipqa[3] = ipqa[2] - 1;
    // printf(" _quad_ck2-ipqa %d %d %d %d / %d %d %d %d\n",
           // ipqa[0],ipqa[1],ipqa[2],ipqa[3], ipqa[0]+MSH2D_pgb_ips,
           // ipqa[1]+MSH2D_pgb_ips,ipqa[2]+MSH2D_pgb_ips,ipqa[3]+MSH2D_pgb_ips);


  return 0;

}


//=========================================================================
  int MSH2D_quad_ckv (int *ifc1, int *en1,
                           int iq1, int iq2, int ip0, int ip1) {
//=========================================================================
// set qma and faces of quads iq1, iq2 to modified and set retCod 
// Input:
//   iq1      quad 
//   iq2      quad right of iq1 / or above iq2
//   qma      1=unmodif; 0=modif.
// Output:
//   ifc1,en1 common-edge
//   retCod   3     both quads unmodified
//            2     iq1 already modified, iq2 unmodified,
//            1     iq1 unmodified, iq2 already modified
//            0     iq1 and iq2 already modified


  int    irc;


  // printf("MSH2D_quad_ckv iq1=%d iq2=%d ip0=%d ip1=%d\n",
          // iq1, iq2, ip0, ip1);
  // printf(" qma[iq1]=%d qma[iq2]=%d\n",MSH2D_qma[iq1],MSH2D_qma[iq2]);

  *ifc1 = -1;

  if(MSH2D_qma[iq1]) {
    // iq1=unmodified: set its faces modified;
    MSH2D_quad_set_mf (ifc1, en1, iq1, ip0, ip1);
    // MSH2D_quad_set_m1 (ifc1, en1, iq1, ip0, ip1);
    MSH2D_qma[iq1] = 0; // set quad = modified
  }


  if(MSH2D_qma[iq2]) {
    // iq2=unmodified: set its faces modified;
    if(*ifc1 < 0) 
      MSH2D_quad_set_mf (ifc1, en1, iq2, ip1, ip0);
    else
      MSH2D_quad_set_m1 (iq2);
    MSH2D_qma[iq2] = 0; // set quad = modified
  }

    // printf("ex-_quad_ckv ifc1=%d en1=%d irc=%d\n",*ifc1,*en1,irc);


  return 0;

}


//=======================================================================
  int MSH2D_quad_set_mf (int *ifc1, int *en1, int iqx, int ip1, int ip2) {
//=======================================================================
// MSH2D_quad_set_m      set quadfaces modified; find edge in quad
// set fac.st of all faces of quad iqx to modified = <MSH2D_qNr>
// Input:
//   iqx        set all faces in this quad modified
//   ip1,ip2    an edge in this quad; find its face#/edge#
// Output:
//   ifc1,en1   face,edge  with points ip1-ip2

  int   i1, fNr, ipe, e1;
  // Fac3  *fa;


  // printf("MSH2D_quad_set_mf iqx=%d ip1=%d ip2=%d\n",iqx,ip1,ip2);

  fNr = MEMTAB_IND(&MSH2D_ts->fac);
  // fa  = MEMTAB_DAT(&MSH2D_ts->fac);

  *ifc1 = -1;

  for(i1=0; i1<fNr; ++i1) {
    if(MFAC__[i1].st != iqx) continue;
    // set face modified
    MFAC__[i1].st = MSH2D_qNr;

    // test if this face has oriented edge ip1-ip2
    // test if this face has point ip1
    e1 = UFA_psn_1fac_ip (&MFAC__[i1], ip1);
    if(e1 < 0) continue;
    ipe = UFA_ipe_fac_esn (&MFAC__[i1], e1);
    if(ipe == ip2) {
      *ifc1 = i1;
      *en1 = e1;
    }
  }

    // printf("ex-_set_mf ifc1=%d en1=%d\n",*ifc1,*en1);

  return 0;

}


//=======================================================================
  int MSH2D_quad_set_m1 (int iqx) {
//=======================================================================
// MSH2D_quad_set_m      set quadfaces modified;
// set fac.st of all faces of quad iqx to modified = <MSH2D_qNr>
// Input:
//   iqx        set all faces in this quad modified
// see also MSH2D_quad_set_mf

  int   i1, fNr;
  // Fac3  *fa;


  // printf("MSH2D_quad_set_m1 iqx=%d\n",iqx);

  fNr = MEMTAB_IND(&MSH2D_ts->fac);
  // fa  = MEMTAB_DAT(&MSH2D_ts->fac);

  for(i1=0; i1<fNr; ++i1) {
    if(MFAC__[i1].st != iqx) continue;
    // set face modified
    MFAC__[i1].st = MSH2D_qNr;
  }

  return 0;

}


//================================================================
  int MSH2D_opt__ () {
//================================================================
// optimize
// was gis_opt__
// 
// was UFA_opt__


  int   irc, i1, fNr, nrf;
  // Fac3  *fa;

  // printf("================================================\n");
  // printf("MSH2D_opt__ start\n");


  L_nxt:

  fNr = MEMTAB_IND(&MSH2D_ts->fac);
  // fa  = MEMTAB_DAT(&MSH2D_ts->fac);

  nrf = 0;


  // loop tru faces
  for(i1=0; i1<fNr; ++i1) {
    // skip unmodified faces (grid-quads, 1) and deleted faces (-1)
    if(MFAC__[i1].st < MSH2D_qNr) continue;
      // printf(" _opt__-f%d\n",i1);


    // optimize face i1, edge 1
    irc = MSH2D_opt_1 (i1, 1);
    if(irc < -90) return irc;
    nrf += irc;

    // optimize face i1, edge 2
    irc = MSH2D_opt_1 (i1, 2);
    if(irc < -90) return irc;
    nrf += irc;

    // optimize face i1, edge 3
    irc = MSH2D_opt_1 (i1, 3);
    if(irc < -90) return irc;
    nrf += irc;

  }
    // printf(" _opt__-nrf = %d\n",nrf);


  if(nrf) goto L_nxt;

  return 0;

}


//================================================================
  int MSH2D_opt_1 (int ifc1, int esn1) {
//================================================================
// optimize iedge esn1 of face ifc1
// returns nr of flips's
// was gis_opt__
// 
//               ip3
//                .    ifc2
//     ip0--------+----------ip2
//                .    ifc1
//               ip1

// was UFA_opt__


  int     irc, i1, fNr, ifc2, esn2, ipa[4], nrf;
  double  d1, d2;
  double  hp1a, hp1b, hp2a, hp2b, hpsa, hpsb;
  Vector2 v01, v02, v03, v12, v13;
  // Fac3    *fna;
  Point2  *pta;


  // printf("MSH2D_opt_1 f%d e%d\n",ifc1,esn1);

  // fa  = MEMTAB_DAT(&MSH2D_ts->fac);
  // fna = MEMTAB_DAT(&MSH2D_ts->fnb);
  pta = MEMTAB_DAT (&MSH2D_ts->pa2);


  //----------------------------------------------------------------
  // get points of contour CCW into ipa

  // get ifc2 = neighbour-face and points
  if(esn1 == 1)      {
    ifc2 = MFNB__[ifc1].i1;
    ipa[0] = MFAC__[ifc1].i1;
    ipa[2] = MFAC__[ifc1].i2;

  } else if(esn1 == 2) {
    ifc2 = MFNB__[ifc1].i2;
    ipa[0] = MFAC__[ifc1].i2;
    ipa[2] = MFAC__[ifc1].i3;

  } else               {
    ifc2 = MFNB__[ifc1].i3;
    ipa[0] = MFAC__[ifc1].i3;
    ipa[2] = MFAC__[ifc1].i1;
  }


  // skip if edge = brkl
  if(ifc2 < 0) return 0;

  // get esn2 = common edge on ifc2
  esn2 = UFA_esn_ips (&MFAC__[ifc2], ipa[2]);
    // printf(" _opt_1 ifc1=%d esn1=%d ifc2=%d esn2=%d ipa[0]=%d ipa[2]=%d \n",
                    // ifc1,esn1, ifc2,esn2, ipa[0],ipa[2]);

  // get ip1 = opposite point of ifc2,esn2
  UFA_ipOpp_fac_esn (&ipa[1], &i1, &MFAC__[ifc2], esn2);

  // get ip3 = opposite point of ifc1,esn1
  UFA_ipOpp_fac_esn (&ipa[3], &i1, &MFAC__[ifc1], esn1);
    // printf(" _opt_1 ipa[1]=%d ipa[3]=%d\n",ipa[1],ipa[3]);


  //----------------------------------------------------------------
  // check convex
  irc = UT2D_is_convex_inpt (4, ipa, pta);
  if(irc >= 0) return 0; // ignore, concave.


  //----------------------------------------------------------------
  // check parameters of p3 and p4 along line p1-p2.
  // parameters must between 0 to 1.
    // printf(" XXXXXXXXXXXX convex: ifc1=%d esn1=%d ifc2=%d esn2=%d \n",
           // ifc1,esn1, ifc2,esn2);
  // get parametric height; if < 0.5 then flip
  // if after flip too acute exit
  // hp1a = parametric-height of ifc1 as is
  // hp1b = parametric-height of ifc1 after flip
  UT2D_vc_2pt (&v13, &pta[ipa[1]], &pta[ipa[3]]);
  UT2D_vc_2pt (&v12, &pta[ipa[1]], &pta[ipa[2]]);
  hp1b = fabs(UT2D_par_nor_2vc (&v13, &v12));
    // printf(" hp1b=%f\n",hp1b);
  if(hp1b < 0.1) return 0;

  UT2D_vc_2pt (&v01, &pta[ipa[0]], &pta[ipa[1]]);
  hp2b = fabs(UT2D_par_nor_2vc (&v13, &v01));
    // printf(" hp2b=%f\n",hp2b);
  if(hp2b < 0.1) return 0;

  UT2D_vc_2pt (&v02, &pta[ipa[0]], &pta[ipa[2]]);
  hp1a = fabs(UT2D_par_nor_2vc (&v02, &v01));

  UT2D_vc_2pt (&v03, &pta[ipa[0]], &pta[ipa[3]]);
  hp2a = fabs(UT2D_par_nor_2vc (&v02, &v03));
    // printf(" _opt_1-flip %d %d\n",ifc1,ifc2);
    // printf("   hp1a=%f hp2a=%f hp1b=%f hp2b=%f\n",hp1a,hp2a,hp1b,hp2b);

  // // if after flip more acute exit
  // if((hp1b + hp2b) < (1.5 * (hp1a + hp2a))) return 0;


  hpsa = DMIN(hp1a, hp2a) + 0.1;
  hpsb = DMIN(hp1b, hp2b);
  if(hpsb <= hpsa) return 0;


  //----------------------------------------------------------------
  // optimize = flip
  // MSH2D_opt_2 ()
  UFA_2fac_flip (ifc1, esn1, ifc2, esn2, MFAC__, MFNB__);

  //----------------------------------------------------------------
  return 1;

}


//================================================================
  int MSH2D_init__ (TessStru *ts1) {
//================================================================
// check if user-defined OB exists; 0=yes, 1=no
// retCod:  >= 0 user-defined OB does exist; index into tab;
//            -1 user-defined OB does NOT exist;
//            -2 ERROR "planar but no OB"

static double fTol = 6.;

  int      irc, i1, itNr, ptNr, fNr, *ia, *ipa, ix1, iy1, ix2, iy2;
  double   dx, dy, d1, d2;
  Point2   p1, p2, *pta;
  IndTab   *itTab;
  GridBox  *qbx;


  // printf("MSH2D_init__ tol=%f\n",ts1->tol);
  // UT3D_stru_dump (Typ_GridBox, &ts1->gbx, ".gbx");


  //----------------------------------------------------------------
  // init
  MSH2D_ts = ts1;
  MFAC__ = MEMTAB_DAT (&MSH2D_ts->fac);
  MFNB__ = MEMTAB_DAT (&MSH2D_ts->fnb);

  // tolerances
  MSH2D_tol = MSH2D_ts->tol;
  MSH2D_2TOL = MSH2D_ts->tol * 2;

  MSH2D_is_delOB = 1;  // 1=faces outside OB exist; 0=not


  // set MSH2D_it_OB >=0: tab-index of OB; -1:  OB-does-not-exist;
  itNr = MEMTAB_IND (&MSH2D_ts->tab);  // nr of contours
  itTab = MEMTAB_DAT (&MSH2D_ts->tab); // contours
  MSH2D_it_OB = -1;
  MSH2D_bnd_nr = 0;
  // loop tru boundaries; find OB
  for(i1=0; i1<itNr; ++i1) {
    if(itTab[i1].typi == MSH_EDGLN_OB) {
      MSH2D_it_OB = i1;
      ++MSH2D_bnd_nr;
    } else if(itTab[i1].typi == MSH_EDGLN_IB) {
      ++MSH2D_bnd_nr;
    }
  }
    // printf(" MSH2D_it_OB=%d MSH2D_bnd_nr=%d\n",MSH2D_it_OB,MSH2D_bnd_nr);


  //----------------------------------------------------------------
  fNr = MEMTAB_IND (&MSH2D_ts->fac);                  // nr of faces
    // printf(" fNr=%d\n",fNr);
  if(fNr > 0) goto L_copy_gbx;


  //================================================================
  // no faces exist yet ..
  // test if OB exists:
  if(MSH2D_it_OB >= 0) {
    // no faces, OB exist;
    // set qbx = gbx reduced (bounding OB)
    MSH2D_init_qbx (&ts1->gbx);
    goto L_exit;

  } else {
    // no faces, no OB
    // check if planar; yes=Error (planar must have OB)
    // if(MSH2D_IS_PLANAR) {
    if(ts1->gbx.ix == 0) {
      TX_Error("MSH2D_tess__ - planar but no OB");
      return -2;
    }
    // copy GridBox
    goto L_copy_gbx;
  }


  //================================================================
  L_copy_gbx:
    // Error if no OB exists
    // qbx = gbx

  // copy GridBox
  MSH2D_qbx = ts1->gbx;

  // set MSH2D_iqxMax = x-index of rightmost quad
  // - this quad is OUTSIDE on right side
  MSH2D_iqxMax = ts1->gbx.ix;

  // set MSH2D_iqyMax = y-index of topmost quad
  // - this quad is OUTSIDE above
  MSH2D_iqyMax = ts1->gbx.iy;


  //----------------------------------------------------------------
  L_exit:

  // check tolerances. dx must be > (tol * 6.)
  d1 = DMIN(MSH2D_qbx.dx, MSH2D_qbx.dy);
  d2 = MSH2D_tol * fTol;
    // printf("     _tess_init %f %f\n",d1,d2);

  if(d1 < d2) {
    MSH2D_tol = d1 / fTol;
    MSH2D_2TOL = MSH2D_ts->tol * 2;
    TX_Print ("****** WNG: tolerance reduced to %f ******",MSH2D_tol);
  }



    // TESTBLOCK
    // printf("MSH2D_iqxMax=%d MSH2D_iqyMax=%d\n",MSH2D_iqxMax,MSH2D_iqyMax);
    // UT3D_stru_dump (Typ_GridBox, &MSH2D_qbx, " qbx-0");
    // MSH2D_test_qbx_disp (qbx);   // display gridlines
    // printf(" ex-_tess_init %d\n",irc);
    // END TESTBLOCK
    
  return irc;

}


//================================================================
  int MSH2D_init_qbx (GridBox *gbx) {
//================================================================
// set qbx = gbx reduced (bounding OB) -
// no faces exist, OB exists.


  int      ptNr, *ipa, *ia, ix1, iy1, ix2, iy2;
  double   dx, dy;
  Point2   p1, p2, *pta;
  IndTab   *itTab;
  GridBox  *qbx;


  // printf("MSH2D_init_qbx \n");

  qbx = &MSH2D_qbx;
  itTab = MEMTAB_DAT (&MSH2D_ts->tab); // contours
  ipa = MEMTAB_DAT (&MSH2D_ts->ipa);
  pta = MEMTAB_DAT (&MSH2D_ts->pa2);   // points

  // find box for OB
  ia = &ipa[itTab[MSH2D_it_OB].ibeg];
  ptNr = itTab[MSH2D_it_OB].iNr;

  // get 2D-bounding-box for pa2-points
  UT2D_box_inpt (&p1, &p2, ptNr, ia, pta);
    // UT3D_stru_dump (Typ_PT2, &p1, "p1");
    // UT3D_stru_dump (Typ_PT2, &p2, "p2");

  // enlarge box 4 * tol 
  UT2D_box_addTol (&p1, &p2, MSH2D_tol * 4);


  //----------------------------------------------------------------
  // planar-surface: ix=0 !  Unlimited support-surface.
  // use 2 faces covering all OB-points.
    if(gbx->ix > 0) goto L_ss1;

    // yes, planar ..
    qbx->ix = 1;
    qbx->iy = 1;

    qbx->pMin.x = p1.x;
    qbx->pMin.y = p1.y;

    qbx->dx = p2.x - p1.x;
    qbx->dy = p2.y - p1.y;

    goto L_exit;


  //----------------------------------------------------------------
  // have dx,dy and ix,iy and pMin from support-surface  (in .gbx)
  L_ss1:

    qbx->dx = gbx->dx;
    qbx->dy = gbx->dy;

    // get ix,iy,quadNr of lower left p1
    DMOD (&ix1, &dx, p1.x, gbx->dx);
    DMOD (&iy1, &dy, p1.y, gbx->dy);
      // printf(" ix1=%d iy1=%d\n",ix1,iy1);

    // get ix,iy,quadNr of upper right p2
    DMOD (&ix2, &dx, p2.x, gbx->dx);
    DMOD (&iy2, &dy, p2.y, gbx->dy);
      // printf(" ix2=%d iy2=%d\n",ix2,iy2);

    // set MSH2D_iqxMax = x-index of rightmost quad
    // - this quad is OUTSIDE on right side
    MSH2D_iqxMax = gbx->ix - ix1;

    // set MSH2D_iqyMax = y-index of topmost quad
    // - this quad is OUTSIDE above
    MSH2D_iqyMax = gbx->iy - iy1;

    // set ix = nr of horiz quads starting at qbx->pMin
    qbx->ix = ix2 - ix1 + 1;
    if((qbx->ix + ix1)  > gbx->ix) qbx->ix = gbx->ix - ix1;

    // set iy = nr of vertic quads starting at qbx->pMin
    qbx->iy = iy2 - iy1 + 1;
    if((qbx->iy + iy1) > gbx->iy) qbx->iy = gbx->iy - iy1;

    qbx->pMin.x = ix1 * qbx->dx;
    qbx->pMin.y = iy1 * qbx->dy;


  //----------------------------------------------------------------
  L_exit:

    // UT3D_stru_dump (Typ_GridBox, &MSH2D_qbx, "ex-init_qbx");

  return 0;

}


//================================================================
  int MSH2D_fac_del_del () {
//================================================================
// MSH_delFac__                 remove deleted faces, update fnb
// delete all faces marked for deletion (fac[].st = 16)
// Overwrite fac[] and fnb[] with last record.
// If (fnb) then update fnb
// was MSH2D_delFac__

  int          irc, i1, ie;
  // Fac3         *fa, *fnba;
  MemTab(Fac3) *facm;


  facm = &MSH2D_ts->fac;
  // fnba = MEMTAB_DAT(&MSH2D_ts->fnb);

  // printf("MSH2D_fac_del_del \n");


  L_pri_ie:
    ie = facm->rNr - 1;    // last record
    if(facm->data[ie].st == -1) {
      facm->rNr -= 1;       // remove last record of facm
      goto L_pri_ie;
    }


  //----------------------------------------------------------------
  // overwrite all deleted records (with valid records)
  for(i1=0; i1<facm->rNr; ++i1) {
    if(facm->data[i1].st != -1) continue;

    // overwrite face[i1] with last face face[ie];
    if(i1 < ie) {
      // update fnb of neighbours
      irc = UFA_fnb_move (ie, i1, facm->data, MFNB__);
      // copy the last fac-record ie -> i1
      facm->data[i1] = facm->data[ie];
      // copy the fnb-record;
      MFNB__[i1] = MFNB__[ie];
        // printf(" _del_del %d -> %d irc=%d\n",ie,i1,irc);   return -99;
    }

    L_nxt_ie:
      facm->rNr -= 1;       // remove last record of facm
      ie = facm->rNr - 1;
      if(facm->data[ie].st == -1) goto L_nxt_ie;
      // --i1;                 // test same face again
      if(ie < 0) { printf("*** MSH2D_fac_del_del E001\n"); return -1; }
  }


/*
  //----------------------------------------------------------------
  // update fnb;
  // fnb-indices are invalid; negative values (breaklines) are correct;
  if(fnba)
  UFA_fnb_if_update (facm->data, fnba, facm->rNr);
*/

  return 0;

}


//================================================================
  int MSH2D_pst_fac_del (Fac3 *fac1) {
//================================================================
// MSH2D_pst_fac_del        mark gridPoints of face i1 unused
// - check all 3 points of face
// - test if point is gridpoint

  int    ipgb;


  // printf("MSH2D_pst_fac_del i1=%d i2=%d i3=%d\n",
         // fac1->i1, fac1->i2, fac1->i3);
  // printf("  ips=%d pNr=%d\n",MSH2D_pgb_ips,MSH2D_mpst.rNr);



  ipgb = fac1->i1 - MSH2D_pgb_ips;    
    // printf(" ipgb-i1=%d\n",ipgb);
  if((ipgb < 0)||(ipgb >= MSH2D_mpst.rNr)) goto L_2;
  if(((char*)MSH2D_mpst.data)[ipgb]) MSH2D_pgst_set (ipgb, 0);


  L_2:
  ipgb = fac1->i2 - MSH2D_pgb_ips;
    // printf(" ipgb-i3=%d\n",ipgb);
  if((ipgb < 0)||(ipgb >= MSH2D_mpst.rNr)) goto L_3;
  if(((char*)MSH2D_mpst.data)[ipgb]) MSH2D_pgst_set (ipgb, 0);


  L_3:
  ipgb = fac1->i3 - MSH2D_pgb_ips;
    // printf(" ipgb-i2=%d\n",ipgb);
  if((ipgb < 0)||(ipgb >= MSH2D_mpst.rNr)) goto L_exit;
  if(((char*)MSH2D_mpst.data)[ipgb]) MSH2D_pgst_set (ipgb, 0);


  L_exit:

  return 0;

}


//================================================================
  int MSH2D_spc_ck1 (int iMore) {
//================================================================
// MSH2D_spc_ck1           check space for additional fac's, fnb's 
//  fnb must have same size as fac


  printf("MSH2D_spc_ck1 %d\n",iMore);

  return 0;

}
 

//================================================================
  int MSH2D_facNrMax (TessStru *ts1) {
//================================================================
// MSH2D_facNrMax               get max nr of faces
// supporting-surface: only space for gridpoints
// bounded-punched-surface: space for gridpoints and breaklines


  int ptNr, fgNr, fbNr;


  // printf("MSH2D_facNrMax \n");


  // nr faces for gridpoints
  fgNr = 2 * ts1->gbx.ix * ts1->gbx.iy;
  fgNr = IMAX (2, fgNr);


  // get nr of used points (boundaries ..)
  fbNr = UFA_fNr_max (MEMTAB_IND (&ts1->pa2));
    // printf(" pa2-rNr,rMax=%d,%d\n",ts1->pa2.rNr,ts1->pa2.rMax);

    // printf("  fgNr=%d fbNr=%d\n",fgNr,fbNr);

  return fgNr + fbNr;

}


//================================================================
  int MSH2D_brkl__ (int ips, int ipe, int edgTyp) {
//================================================================
// make breakline of boundary-segment ips-ipe
// Both endpoints of breakline are already meshed.
// Input:
//   edgTyp         -4 = MSH_EDGLN_OB | -3 = MSH_EDGLN_IB
// retCod           0 OK;
//                 -1 ERR
// was MSH_msh_BL__


  int      irc, i1, i2, i3, i4, fNr, fnb1, fnb2, esn1, esn2, siz_cfa, c1Nr, c2Nr;
  int      *cfa, *c1fa, *c2fa, *dfa, *fan, nf1a, ifc1, ifc2;
  int      ecNr, fiNr, *fia;
  int      *ipba, *foa, *eoa;
  char     *cea, *c1ea, *c2ea;
  Fac3     *fc1a, *fc2a;
  // ecStru   *eca;
  brklStru brkls;
  // Point2   *pta;
  MemTab(Fac3) f1a=_MEMTAB_NUL;
  TessStru *ts1 = MSH2D_ts;


  // printf("IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII \n");
  // printf("MSH2D_brkl__ ips=%d ipe=%d edgTyp=%d\n",ips,ipe,edgTyp);
  // TESTBLOCK
  // if(ips == 1) return -99;
  // END TESTBLOCK


  // fa  = MEMTAB_DAT(&ts1->fac);
  fNr = MEMTAB_IND(&ts1->fac);
  // fnb = MEMTAB_DAT(&ts1->fnb);
  // pta = MEMTAB_DAT (&MSH2D_ts->pa2);

  // get space for list of edges (cfa,cea)
  siz_cfa = MSH2D_edgnr_facnr(fNr) / 2;  
    // printf(" siz_cfa=%d\n",siz_cfa);


  //----------------------------------------------------------------
  brkls.c1Nr  = 0;
  brkls.c2Nr  = 0;
  brkls.fiNr   = 0;
  brkls.ips    = ips;
  brkls.ipe    = ipe;
  brkls.siz    = siz_cfa;
  brkls.edgTyp = edgTyp;
  brkls.pta    = MEMTAB_DAT (&MSH2D_ts->pa2);
  brkls.ps     = &brkls.pta[ips];
  brkls.pe     = &brkls.pta[ipe];
    // printf(" ps = %f %f\n",brkls.ps->x,brkls.ps->y);
    // printf(" pe = %f %f\n",brkls.pe->x,brkls.pe->y);
  UT2D_vc_2pt (&brkls.vcl, brkls.ps,  brkls.pe);

  // get space for right contour
  brkls.cf1a = (int*) MEM_alloc_tmp (siz_cfa * (int)sizeof(int));
  if(brkls.cf1a == NULL) {TX_Error("MSH2D_brkl__ EOM-1"); return -1;}

  // get space for left contour
  brkls.cf2a = (int*) MEM_alloc_tmp (siz_cfa * (int)sizeof(int));
  if(brkls.cf2a == NULL) {TX_Error("MSH2D_brkl__ EOM-2"); return -1;}

  // get space for list of faces to delete (inside contour)
  brkls.fia = (int*) MEM_alloc_tmp (siz_cfa * (int)sizeof(int));
  if(brkls.fia == NULL) {TX_Error("MSH2D_brkl__ EOM-3"); return -1;}

  // get space for edgNr for cf1a
  brkls.ce1a = (char*) MEM_alloc_tmp (siz_cfa);
  if(brkls.ce1a == NULL) {TX_Error("MSH2D_brkl__ EOM-4"); return -1;}

  // get space for edgNr for cf1a
  brkls.ce2a = (char*) MEM_alloc_tmp (siz_cfa);
  if(brkls.ce2a == NULL) {TX_Error("MSH2D_brkl__ EOM-5"); return -1;}


  //----------------------------------------------------------------
  // get bounding-contours right and left of brkLn ips-ipe;
  //   and add all faces to inner-faces-list
  irc = MSH2D_brkl_bnd__ (&brkls);
  if(irc < -90) return irc;

    // TESTBLOCK
    // if(brkls.ips == 12) {
      // MSH2D_dump_brkls (&brkls, "b");
      // UFA_Disp_nEdg2 (brkls.fiNr, brkls.fia, NULL,
                      // MFAC__, brkls.pta, "s", 0, -0.001);
      // return -99;
    // }
    // END TESTBLOCK


  //----------------------------------------------------------------
  if(irc == 1) {
    // breakline already exists
    MSH2D_brkl_set (brkls.ixf, brkls.ixe, ips, ipe, edgTyp);
    goto L_ffin;
  }

  //----------------------------------------------------------------
  // only one crossing-edge; flip ?
  // if(brkls.c1Nr == 1) {
  if(irc == 2) {
    irc = MSH2D_brkl_flip (brkls.ixf, brkls.ixe, ips, ipe, edgTyp);
    goto L_ffin;
  }


  //----------------------------------------------------------------
  // get space for ipba = list of boundaryPoints;
  i1 = brkls.c1Nr + brkls.c2Nr + 4;
  ipba = (int*) MEM_alloc_tmp (i1 * (int)sizeof(int));
  if(ipba == NULL) {TX_Error("MSH2D_brkl__ EOM-3A"); return -1;}

  // get space for new faces
  // i1 = MSH2D_facnr_edgnr(c1Nr) + MSH2D_facnr_edgnr(c2Nr) + 4;
  i1 = brkls.c1Nr + brkls.c2Nr + 4;
  i1 *= sizeof(Fac3);
  fc1a = (Fac3*) MEM_alloc_tmp (i1);
  if(fc1a == NULL) {TX_Error("MSH2D_brkl__ EOM-31"); return -1;}
  MemTab_spc (&f1a, fc1a, i1, sizeof(Fac3), Typ_Fac3);

  // reverse brkls.cf2a, brkls.ce2a
  if(brkls.c2Nr > 0) {
    MEM_inv_itab (brkls.c2Nr, brkls.cf2a);
    MEM_inv_btab (brkls.c2Nr, brkls.ce2a);
  }
    // for(i1=0;i1<brkls.c2Nr;++i1)
      // printf(" c2[%d] f%d e%d\n",i1,brkls.cf2a[i1],brkls.ce2a[i1]);


  // create ipba = list of boundaryPoints
  irc = MSH2D_brkl_ipba (ipba, &brkls);
  if(irc < -90) return irc;

    // TESTBLOCK
    // if(ips == 6) {
      // MSH2D_view_cv_inpt (brkls.c1Nr + brkls.c2Nr + 1, ipba);
      // return -99;
    // }
    // END TESTBLOCK


  // get i1 = nr of points for right contour
  i1 = brkls.c1Nr + 1;

  // get i2 = nr of points for left contour
  i2 = brkls.c2Nr + 1;

  // tesselate ips-ipe
  // irc = MSH2D_ctess_e (&f1a, brkls.c1Nr, brkls.cf1a, brkls.ce1a);
  if(i1 > 1) {
    irc = MSH2D_ctess__ (&f1a, ipba, i1, brkls.pta);
      // UFA_fac_dump_ (MEMTAB_DAT(&f1a), MEMTAB_IND(&f1a), " f1a-1");
    if(irc < -90) return irc;
  }

  // tesselate ipe-ips
  // irc = MSH2D_ctess_e (&f1a, brkls.c2Nr, brkls.cf2a, brkls.ce2a);
  if(i2 > 1) {
    irc = MSH2D_ctess__ (&f1a, &ipba[brkls.c1Nr], i2, brkls.pta);
      // UFA_fac_dump_ (MEMTAB_DAT(&f1a), MEMTAB_IND(&f1a), " f1a-2");
    if(irc < -90) return irc;
  }

    // TESTBLOCK
    // if(ips == 0) {
      // UFA_view2__ (MEMTAB_DAT(&f1a), 0, MEMTAB_IND(&f1a),
                   // brkls.pta, MEMTAB_IND (&MSH2D_ts->pa2), 
                   // "sfb", -1L, 0.1);
      // return -99;
    // }
    // END TESTBLOCK



  // set all quads of faces to delete to modified !
  if(MSH2D_qma) {
    irc = MSH2D_quad_nf_mod (brkls.fiNr, brkls.fia);
    if(irc < -90) return irc;
  }



  //----------------------------------------------------------------
  // get memSpc for outer-neighbours-of-boundary
  i1 = brkls.c1Nr + brkls.c2Nr;
  foa = (int*) MEM_alloc_tmp (i1 * (int)sizeof(int));
  eoa = (int*) MEM_alloc_tmp (i1 * (int)sizeof(int));

  // get outer-neighbours-of-boundary
  irc = MSH2D_brkl_fnb_eo (foa, eoa, &brkls);
  if(irc < 0) return irc;


    // TESTBLOCK
    // if(brkls.ips == 6) {
      // UFA_Disp_nEdg2 (i1,foa,eoa,MFAC__,brkls.pta,"bs",Typ_Att_hili1,0.1);
      // return -99;
    // }
    // END TESTBLOCK



  //----------------------------------------------------------------
  // delete all faces inside contour ips-ipe-ips
  // irc = MSH2D_fdel_bnd (c1Nr + c2Nr, cfa, cea);
  // if(irc < -90) return irc;
  for(i1=0; i1 < brkls.fiNr; ++i1) MFAC__[brkls.fia[i1]].st = -1;


  //----------------------------------------------------------------
  // get space for list of new faces
  nf1a = MEMTAB_IND (&f1a);
  fan = (int*) MEM_alloc_tmp (nf1a * (int)sizeof(int));
    // printf(" nf1a=%d\n",nf1a);


  // add f1a = new faces to ts1.fac (overwriting deleted faces)
  // get fan = list of new faces
  // remove all deleted faces
  irc = MSH2D_brkl_cpy (fan, &f1a);
  if(irc < -90) return irc;


  // update fnb of new-faces and of neighbours of new-faces
  irc = MSH2D_brkl_fnb__ (nf1a, fan, foa, eoa, ipba, &brkls);
  if(irc < 0) return irc;


  //----------------------------------------------------------------
  L_ffin:
    // TESTBLOCK
    // printf("ex MSH2D_brkl__ ips=%d ipe=%d edgTyp=%d\n",ips,ipe,edgTyp);
    // MSH2D_view_nfac ("bf", 0, -1, "OB_1");
    // irc = MSH2D_dump__ ("ex-brkl__"); if(irc) return irc;
    // if(ips == 6) {
      // irc = UFA_ifev_find_2ip (&i1,&i2,4,3,MEMTAB_DAT(&ts1->fac),
                               // MSH2D_fnb,MEMTAB_IND(&ts1->fac));
      // return -99;
    // }
    // END TESTBLOCK

    return 0;

}


//================================================================
  int MSH2D_brkl_ipba (int *ipba, brklStru *brkls) {
//================================================================
// create ipba = list of boundaryPoints (CCW, around right and left contour)
//  c1: <brkls->c1Nr + 1> points, first = 
//  c2: <brkls->c2Nr + 1> points, first = <brkls->c1Nr>
// was MSH2D_ctess_e

  int   i1, i2, ii, *cfa;
  char  *cea;


  // printf("MSH2D_brkl_ipba ips=%d\n",brkls->ips);

  ii = 0;


  // contour-right
  if(brkls->c1Nr > 0) {
    // add startpoint of all edges of contour-right
    i2 = brkls->c1Nr;
    cfa = brkls->cf1a;
    cea = brkls->ce1a;
    for(i1=0; i1<i2; ++i1) {
      ipba[ii] = UFA_ips_fac_esn (&MFAC__[cfa[i1]], cea[i1]);
      ++ii;
    }
  }

  // contour-left
  if(brkls->c2Nr > 0) {
    // add startpoint of all edges of contour-left
    i2 = brkls->c2Nr;
    cfa = brkls->cf2a;
    cea = brkls->ce2a;
    for(i1=0; i1<i2; ++i1) {
      ipba[ii] = UFA_ips_fac_esn (&MFAC__[cfa[i1]], cea[i1]);
      ++ii;
    }
  }

  // add endpoint of last edge
  --i2;
  ipba[ii] = UFA_ipe_fac_esn (&MFAC__[cfa[i2]], cea[i2]);
  ++ii;

  // add startpoint of first edge
  ipba[ii] = ipba[0];


    // TESTBLOCK
    // printf("ex-MSH2D_brkl_ipba %d\n",brkls->c1Nr);
    // i2 = brkls->c1Nr + brkls->c2Nr + 1;
    // for(i1=0;i1<i2;++i1) printf(" ipba[%d] = %d\n",i1,ipba[i1]);
    // printf(" c1: %d points, first = 0\n",brkls->c1Nr + 1);
    // printf(" c2: %d points, first = %d\n", brkls->c2Nr + 1, brkls->c1Nr);
    // if(brkls->ips == 0) {
      // MSH2D_view_cv_inpt (brkls->c1Nr + brkls->c2Nr + 1, ipba);
      // return -99;
    // }
    // END TESTBLOCK



  return 0;

}


//================================================================
  int MSH2D_brkl_fnb_eo (int *foa, int *eoa, brklStru *brkls) {
//================================================================
// get outer-neighbours-of-boundaries in brkls
// Output: foa,eoa = edges of neighbours-faces around cf1a,cf2a-boundary

  
  int      irc, i1, i2, ii;


  // printf("MSH2D_brkl_fnb_eo ----------\n");
  // MSH2D_dump_brkls (brkls, "");

  // get outer-neighbours-of-boundary
  ii = 0;
  i2 = brkls->c1Nr;
  for(i1=0; i1<i2; ++i1) {
    irc = UFA_ife_getNf_ife (&foa[ii], &eoa[ii],
                             brkls->cf1a[i1], brkls->ce1a[i1], MFAC__, MFNB__);
    ++ii;
  }

  i2 = brkls->c2Nr;
  for(i1=0; i1<i2; ++i1) {
    irc = UFA_ife_getNf_ife (&foa[ii], &eoa[ii],
                             brkls->cf2a[i1], brkls->ce2a[i1], MFAC__, MFNB__);
    ++ii;
  }

    // TESTBLOCK
    // i2 = brkls->c1Nr + brkls->c2Nr;
    // for(i1=0;i1<i2;++i1) printf(" _fnb_eo-fo[%d] %d %d\n",i1,foa[i1],eoa[i1]);
    // if(brkls->ips == 4)
      // UFA_Disp_nEdg2 (i2,foa,eoa, MFAC__,brkls->pta, "bs",Typ_Att_hili1,0.);
    // if(brkls->ips == 4)
      // return -99;
    // END TESTBLOCK


  return 0;

}


//================================================================================
  int MSH2D_brkl_fnb__ (int fiNr, int *fia, int *foa, int *eoa,
                        int *ipba, brklStru *brkls) {
//================================================================================
// update fnb of new-faces and of neighbours of new-faces

  int      irc, i1, i2, ii, foNr, ips, ipe, if1, en1;
  // Fac3     *fnb;


  // printf("MSH2D_brkl_fnb__ %d\n",fiNr);
  // MSH2D_dump_brkls (brkls, "");
  // for(i1=0; i1<fiNr; ++i1) printf(" fia[%d] = %d\n",i1,fia[i1]); // inner faces
  // i2 = brkls->c1Nr + brkls->c2Nr + 1;
  // for(i1=0;i1<i2;++i1) printf(" ipba[%d] = %d\n",i1,ipba[i1]);


  // fa = MEMTAB_DAT (&MSH2D_ts->fac);
  // fnb = MEMTAB_DAT(&MSH2D_ts->fnb);


  //----------------------------------------------------------------
  // loop tru outer-boundary; set fnb
  foNr = brkls->c1Nr + brkls->c2Nr;
  for(i1=0; i1<foNr; ++i1) {
      ips = ipba[i1];
      ipe = ipba[i1+1];
      // printf(" nxt-outer[%d] f%d e%d  ips=%d ipe=%d\n",
             // i1,foa[i1],eoa[i1],ips,ipe);
    // check for breakline
    if(foa[i1] < 0) {
      // find (inner) edge with ips-ipe in fia
      irc = UFA_if_nf_2ip (&if1, &en1, fiNr, fia, ips, ipe, MFAC__);
      if(irc < 0) return irc;
      // set fnb on outer face to breakline
      UFA_fnb_set_if_es (MFNB__, if1, en1, foa[i1]);
      continue; // outer-edge foa[i1],eoa[i1] is a breakline
    }


    // find inner-edge ipe-ips in fia (inner faces)
    irc = UFA_if_nf_2ip (&if1, &en1, fiNr, fia, ips, ipe, MFAC__);
    if(irc < 0) return irc;
      // printf(" innerEdge[%d]=f%d,e%d irc=%d\n",i1,if1,en1,irc);

    // set fnb on outer face
    UFA_fnb_set_if_es (MFNB__, foa[i1], eoa[i1], if1);

    // set fnb on inner face
    UFA_fnb_set_if_es (MFNB__, if1, en1, foa[i1]);
  }


  //----------------------------------------------------------------
  // set breakline on inner faces
    // printf(" ---------------brkl_fnb__-2\n");
    // printf(" ips=%d ipe=%d\n",brkls->ips,brkls->ipe);
  // find edge with ips-ipe in fia
  irc = UFA_if_nf_2ip (&if1, &en1, fiNr, fia, brkls->ips, brkls->ipe, MFAC__);
  // if breakline covers outer boundary, no edge might exist
  // if(irc < 0) return irc;
  if(!irc) UFA_fnb_set_if_es (MFNB__, if1, en1, brkls->edgTyp);

  // find edge with ipe-ips in fia
  irc = UFA_if_nf_2ip (&if1, &en1, fiNr, fia, brkls->ipe, brkls->ips, MFAC__);
  // if(irc < 0) return irc;
  if(!irc) UFA_fnb_set_if_es (MFNB__, if1, en1, brkls->edgTyp);




  //----------------------------------------------------------------
  // set (missing) fnb for all inner faces
  // loop tru all inner faces
    // printf(" ---------------brkl_fnb__-3\n");

  for(i1=0; i1<fiNr; ++i1) {
    ii = fia[i1];
      // UFA_fnb_dump_1 (&MFAC__[ii], &MFNB__[ii], "inner[%d] = f%d ",i1,ii);

    if(MFNB__[ii].i1 == MSH_EDGLN_ERR) {
      ips = MFAC__[ii].i2;
      ipe = MFAC__[ii].i1;
      // find (inner) edge with ips-ipe in fia
      irc = UFA_if_nf_2ip (&if1, &en1, fiNr, fia, ips, ipe, MFAC__);
      if(irc < 0) return irc;
      // set fnb on inner face
      UFA_fnb_set_if_es (MFNB__, ii, 1, if1);
      UFA_fnb_set_if_es (MFNB__, if1, en1, ii);
    }
        
    if(MFNB__[ii].i2 == MSH_EDGLN_ERR) { 
      ips = MFAC__[ii].i3;
      ipe = MFAC__[ii].i2;
      // find (inner) edge with ips-ipe in fia
      irc = UFA_if_nf_2ip (&if1, &en1, fiNr, fia, ips, ipe, MFAC__);
      if(irc < 0) return irc;
      // set fnb on inner face
      UFA_fnb_set_if_es (MFNB__, ii, 2, if1);
      UFA_fnb_set_if_es (MFNB__, if1, en1, ii);
    }

    if(MFNB__[ii].i3 == MSH_EDGLN_ERR) { 
      ips = MFAC__[ii].i1;
      ipe = MFAC__[ii].i3;
      // find (inner) edge with ips-ipe in fia
      irc = UFA_if_nf_2ip (&if1, &en1, fiNr, fia, ips, ipe, MFAC__);
      if(irc < 0) return irc;
      // set fnb on inner face
      UFA_fnb_set_if_es (MFNB__, ii, 3, if1);
      UFA_fnb_set_if_es (MFNB__, if1, en1, ii);
    }

    // set st=OK
    MFNB__[ii].st = 0;

  }
  


  return 0;

}

 
//===================================================================
  int MSH2D_brkl_bnd__ (brklStru *brkls) {
//===================================================================
// find all edges crossing line ips-ipe
// Input:
//   eCxNr      size of eCxa
//   ips,ipe    line
// Output:
//   eCxa       table of crossing edges
//   RetCod:
//      0       OK; contours & faces in brkls
//      1       breakline already exists; ifc in 1. and 2. eCxa-record
//      2       flip brkls->ixf,ixe
//
// was MSH2D_nife_segCov


  int      irc, i1;
  // Point2   *ps, *pe, *pta;


    // TESTBLOCK
    // if(brkls->ipe == 2) return -99;
    // END TESTBLOCK


  // fNr = MEMTAB_IND (&MSH2D_ts->fac);
  // fa = MEMTAB_DAT (&MSH2D_ts->fac);

  // pta = MEMTAB_DAT (&MSH2D_ts->pa2);
  // ps = &pta[ips];
  // pe = &pta[ipe];
    // printf(" ps = %f %f\n",ps->x,ps->y);
    // printf(" pe = %f %f\n",pe->x,pe->y);


  // get first intersecting edge
  irc = MSH2D_brkl_bnd_1_ (brkls);
  if(irc != 0) return irc;


  //----------------------------------------------------------------
  // add all following intersecting edges
  L_nxt_0:
  irc = MSH2D_brkl_bnd_2_ (brkls);
  if(irc < 0) return irc;



  //----------------------------------------------------------------
  L_exit:
    // TESTBLOCK
    // printf("ex-MSH2D_brkl_bnd__ irc=%d\n",irc);
    // if(brkls->ipe == 4) {MSH2D_dump_brkls (brkls, "bs");return -99;}
    // MSH2D_dump_brkls (brkls, "ex-MSH2D_brkl_bnd__");
    // END TESTBLOCK


  return irc;

}


//=====================================================================
  int MSH2D_brkl_bnd_11 (brklStru *brkls, int if1, int en1, int mode) {
//=====================================================================
// edge if1,en1 is parallel ips-ipe;
// Input:
//   mode    0 if1,en1 is right of ips-ipe; ends at line
//           1 if1,en1 is left of ips-ipe; starts at line.


  int     irc, enp1, enn1, if2, en2, enp2, enn2;


  // printf("MSH2D_brkl_bnd_11 f%d, e%d mode=%d\n",if1,en1,mode);


  //----------------------------------------------------------------
  if(mode == 0) {
    // if1,en1 is right of ips-ipe; ends at ips-ipe

    // get prv & nxt for if1
    UFA_2esn_prvnxt (&enp1, &enn1, en1);

    // get if2 = nbf of enn1
    UFA_ife_getNf_ife (&if2, &enp2, if1, enn1, MFAC__, MFNB__);
      // printf(" m0: if2=%d enp2=%d\n",if2,enp2);

    // is edge on outer boundary ?
    if(enp2 < 0) goto L_right_only;


    // get prv & nxt for if2
    UFA_2esn_prvnxt (&enn2, &en2, enp2);


  //----------------------------------------------------------------
  } else {
    // if1,en1 is left of ips-ipe; starts at ips-ipe
    if2 = if1;
    en2 = en1;

    // get prv & nxt for if2
    UFA_2esn_prvnxt (&enp2, &enn2, en2);

    // get if2 = nbf of edge ips-iee
    UFA_ife_getNf_ife (&if1, &enn1, if2, enp2, MFAC__, MFNB__);
      // printf(" m1: if1=%d enn1=%d\n",if2,enn1);

    // is edge on outer boundary ?
    if(enn1 < 0) goto L_left_only;

    // get prv & nxt for if1
    UFA_2esn_prvnxt (&en1, &enp1, enn1);
  }


  //----------------------------------------------------------------
    // printf(" en1=%d enn1=%d enp1=%d\n",en1,enn1,enp1);
    // printf(" en2=%d enn2=%d enp2=%d\n",en2,enn2,enp2);


  // add both faces to list
  // add if1 to facelist
  irc = MSH2D_brkl_bnd_addf (brkls, if1);
  if(irc < 0) return -1;

  // add if2 to facelist
  irc = MSH2D_brkl_bnd_addf (brkls, if2);
  if(irc < 0) return -1;

  // add previous edge of if1,en1 to right-boundary and set brkls->cp1e
  irc = MSH2D_brkl_bnd_add1 (brkls, if1, enp1);
  if(irc < 0) return irc;

  // and add next edge to left-boundary and set brkls->cp2s
  irc = MSH2D_brkl_bnd_add2 (brkls, if2, enn2);
  if(irc < 0) return irc;

  // set ixf,ixe = edge from end of contour-right (cp1e) to point on ips-ipe
  brkls->ixf  = if1;
  brkls->ixe  = en1;   // last edge right of ips-ipe (cp1e-ipxf)
  brkls->ixlf = if2;
  brkls->ixle = en2;   // last edge left of ips-ipe (cp2s-ipxf)
  brkls->ixTyp = 1;    // point is on edge ips-ipe


  return 0;


  //================================================================
  L_right_only:
  // edge if1,enn1 is on/parallel ips-ipe;
  brkls->ixTyp = 3;    // point is on edge ips-ipe
  brkls->ixf  = if1;
  brkls->ixe  = en1;   // last edge right of ips-ipe (cp1e-ipxf)

  // add if1 to facelist
  irc = MSH2D_brkl_bnd_addf (brkls, if1);
  if(irc < 0) return -1;

  // add previous edge of if1,en1 to right-boundary and set brkls->cp1e
  irc = MSH2D_brkl_bnd_add1 (brkls, if1, enp1);
  if(irc < 0) return irc;

  // // add next edge of if1,en1 to left-boundary and set brkls->cp2s
  // irc = MSH2D_brkl_bnd_add2 (brkls, if1, enn1);
  // if(irc < 0) return irc;

  return 0;


  //================================================================
  L_left_only:
  // edge if2,enp2 is on/parallel ips-ipe;
  brkls->ixTyp = 4;    // point is on edge ips-ipe
  brkls->ixlf = if2;
  brkls->ixle = en2;   // last edge left of ips-ipe (cp2s-ipxf)

  // add if2 to facelist
  irc = MSH2D_brkl_bnd_addf (brkls, if2);
  if(irc < 0) return -1;

  // // add prev edge to right-boundary and set brkls->cp1e
  // irc = MSH2D_brkl_bnd_add1 (brkls, if2, enp2);
  // if(irc < 0) return irc;

  // add next edge to left-boundary and set brkls->cp2s
  irc = MSH2D_brkl_bnd_add2 (brkls, if2, enn2);
  if(irc < 0) return irc;

  return 0;

}


//===================================================================
  int MSH2D_brkl_bnd_1_ (brklStru *brkls) {
//===================================================================
// find first edge crossing line ips-ipe
// Input:
//   brkls->ips,ipe  line
// Output:
//   brkls->cf1a      1. (edge) right of ips-ipe
//   brkls->cf2a      1. (edge) left of ips-ipe
//   brkls->fia       1. face intersecting ips-ipe
//   brkls->ixTyp     0 = normal intersection;
//                    1 = edge is parallel ips-ipe; face = right
//                    2 = edge is parallel ips-ipe; face is left
//   brkls->ixf,ixe   ixTyp-0: crossing edge cp1e-cp2s
//                    ixTyp-1: last edge cp1e-ipxf
//   brkls->ixlf,ixle ixTyp-1: last edge ipxf-cp2s
//   brkls->cp1e      endPt of right-contour
//   brkls->cp2s      startPt of left-contour
//   brkls->ipxf      ixTyp-1: point of ixf/ixlf on ips-ipe
//   RetCod:
//      0       1. crossing-edges in brkls;
//      1       breakline already exists; ifc in 1. eca-record
//      2       only one crossing-edge; flip
//
// was MSH2D_nife_segCov


#define PSTABSIZ 100

  int      irc, iNr, i1, lTyp, fNr, if1, en1, ies, iee, if2, en2, ipo;
  int      enp1, enn1, enp2, enn2;
  int      iTab[PSTABSIZ];
  char     psna[PSTABSIZ];
  Point2   *p21, *p22, *pta;
  // Fac3     *fnb;


  // printf("MSH2D_brkl_bnd_1_ ips=%d ipe=%d\n",brkls->ips,brkls->ipe);
  // if(brkls->ips == 2) return -99;


  fNr = MEMTAB_IND (&MSH2D_ts->fac);
  // fa = MEMTAB_DAT (&MSH2D_ts->fac);
  // fnb = MEMTAB_DAT(&MSH2D_ts->fnb);

  pta = brkls->pta;


  //----------------------------------------------------------------
  // get iTab = table of all faces using point ips.
  // And get psna = edge of face iTab starting with point ips 
  iNr = PSTABSIZ;
  irc = UFA_nifac_ck_pt (iTab, psna, &iNr, brkls->ips, MFAC__, fNr);
  if(irc < 0) return -1;

    // TESTBLOCK
    // for(i1=0;i1<iNr;++i1)
      // UFA_fac_dump1 (&MFAC__[iTab[i1]], "%d-e%d-touch-ips",iTab[i1],psna[i1]);
    // END TESTBLOCK


  //----------------------------------------------------------------
  // test if edge ips-ipe already exists.
  for(i1=0; i1<iNr; ++i1) {
    // test if endpoint of edge == ipe
    if(brkls->ipe == UFA_ipe_fac_esn (&MFAC__[iTab[i1]], psna[i1])) {
      // yes, irc=1
      brkls->ixf = iTab[i1];
      brkls->ixe = psna[i1];
      irc = 1;
      goto L_exit;
    }
  }


  //----------------------------------------------------------------
  // find first edge intersecting line
  // test intersection fc1, es1 with ips-ipe
  for(i1=0;i1<iNr;++i1) {
    if1 = iTab[i1];
    en1 = UFA_esn_nxt (psna[i1]);

    // get start,endPoint
    UFA_ipse_fac_esn (&ies, &iee, &MFAC__[if1], en1);
    p21 = &pta[ies];
    p22 = &pta[iee];
      // printf(" if1=%d en1=%d ies=%d iee=%d\n",if1,en1,ies,iee);

    // intersect brlLn (ps-pe) with edge p21-p22
    irc = UT2D_2pt_ck_int_2pt (&lTyp, &MSH2D_tol, brkls->ps, brkls->pe, p21, p22);
      // printf("   ex-ck_int_2pt-i1=%d irc=%d lTyp=%d\n",i1,irc,lTyp);
    if(irc < 0) continue;   // no intersection ..


    //----------------------------------------------------------------
    if(irc < 1) {
      // 0 = normal intersection;
      brkls->ixTyp = 0;
      brkls->ixf = if1;
      brkls->ixe = en1;
      // get prev & nxt
      UFA_2esn_prvnxt (&enp1, &enn1, en1);
      // add previous edge to right-boundary and set brkls->cp1e
      irc = MSH2D_brkl_bnd_add1 (brkls, if1, enp1);
      if(irc < 0) return irc;
      // and add next edge to left-boundary and set brkls->cp2s
      irc = MSH2D_brkl_bnd_add2 (brkls, if1, enn1);
      if(irc < 0) return irc;
      // and add if1 to facelist
      irc = MSH2D_brkl_bnd_addf (brkls, if1);
      if(irc < 0) return -1;
      goto L_exit;
    }


    //----------------------------------------------------------------
    if(irc == 1) {
      // 1   edge starts at line ips-ipe; face is left
      // mark point ies unused
      MSH2D_pst_set (ies, 0);
      // point of face ixf/ixlf on ips-ipe
      brkls->ipxf = ies;
      irc = MSH2D_brkl_bnd_11 (brkls, if1, en1, 1);
      if(irc < 0) return irc;
      goto L_exit;
    }


    //----------------------------------------------------------------
    if(irc == 2) {
      // 2   endPt iee=p22 of if1,en1 is on line ps-pe; face = right
      // mark point iee unused
      MSH2D_pst_set (iee, 0);
      // point of face ixf/ixlf on ips-ipe
      brkls->ipxf = iee;
      irc = MSH2D_brkl_bnd_11 (brkls, if1, en1, 0);
      if(irc < 0) return irc;
      goto L_exit;
    }
 

    //----------------------------------------------------------------
      // 3 (p1=p3), 4 (p1=p4), 5 (p2=p3), 6 (p2=p4) identical endpoints
      // 7 = parallel lines with overlapping

    if(irc == 5) {
      // p2=p3; ipe=ies; 
    }

TX_Error("MSH2D_brkl_bnd_1_ E-1-%d",irc); return -99;

  }



  //----------------------------------------------------------------
  L_exit:
    // TESTBLOCK
    // printf("ex-MSH2D_brkl_bnd_1_ irc=%d\n",irc);
    // MSH2D_dump_brkls (brkls, "ex-MSH2D_brkl_bnd");
    // if(brkls->ips == 0)
    // return -99;
    // END TESTBLOCK


  return irc;

}


//===================================================================
  int MSH2D_brkl_bnd_2_ (brklStru *brkls) {
//===================================================================
// find all edges crossing line ips-ipe
// add contours
// add faces 
// Input:
//   brkls
// Output:
//   RetCod:
//      0       OK; contours & faces in brkls
//      2       flip brkls->ixf,ixe
//      9       unfinished
//
// was MSH2D_nife_segCov


  int      irc, i1, if1, en1, if2, en2, spe, ixTyp, enp2, enn2;
  int      ixps, ixpe, ip2s, ip2e, ip2o, ipll, ipxf;
  int      ixlf,ixle, if3, en3, enp3, ixNxt;
  Point2   *pta;


  // printf("MSH2D_brkl_bnd_2_ ips=%d ipe=%d\n",brkls->ips,brkls->ipe);
  // if(brkls->ips == 0) return -99;

  pta = brkls->pta;

  ixTyp = brkls->ixTyp;



  // last crossing-edge
  if1 = brkls->ixf;
  en1 = brkls->ixe;

  ixlf = brkls->ixlf;
  ixle = brkls->ixle;

  // pt on ips-ipe; endPt of if1,en1; startPt of ixlf,ixle (only for ixTyp=1)
  ipxf = brkls->ipxf;
  

  //----------------------------------------------------------------
  L_nxt:
    // printf("L_nxt: -------------------------------------- \n");

      // TESTBLOCK
      // printf(" ixTyp=%d if1=%d en1=%d cp1e=%d cp2s=%d\n",
              // ixTyp,if1,en1,brkls->cp1e,brkls->cp2s);
      // if(if1 == 28) return -99;
      // END TESTBLOCK


  //================================================================
  if(brkls->ixTyp == 3) {
    // startPt ipxf is on line ips-ipe;
    // ixlf,ixle is ipxf-cp1e
    irc = MSH2D_brkl_bnd_23 (brkls);
    goto L_exit;


  //================================================================
  } else if(brkls->ixTyp == 4) {
    // startPt ipxf is on line ips-ipe;
    // ixlf,ixle is ipxf-cp2s
    irc = MSH2D_brkl_bnd_24 (brkls);
    goto L_exit;
  }


  //================================================================
  if(ixTyp != 0) goto L_typ_1_;
  // last was normal crossing-edge
  // Using: if1,en1 ist the edge crossing the line;
 
  // get neigbour-face and common edge of if1,en1
  UFA_ife_getNf_ife (&if2, &enn2, if1, en1, MFAC__, MFNB__);
  if(if2 < 0) {TX_Error("MSH2D_brkl_bnd_2_-E2");return -91;}

  // get enp2, enn2
  UFA_2esn_prvnxt (&en2, &enp2, enn2);

  // add face
  irc = MSH2D_brkl_bnd_addf (brkls, if2);
  if(irc < 0) return -1;

  // get ip2o = the opposite-point of if2,en2
  ip2o = UFA_ipo_fac_esn (&MFAC__[if2], enn2);
    // printf("  L_typ_00: if2=%d en2=%d ip2o=%d\n",if2,en2,ip2o);

  // test if ip2o == ipe
  if(ip2o == brkls->ipe) {
    // opposite-point of last crossing-edge = endpoint; done.
    // if only one normal intersection: special-case-flip.
    if(brkls->fiNr == 1) { irc = 2; goto L_exit; }

    // add prev edge as right-contour
    irc = MSH2D_brkl_bnd_add1 (brkls, if2, enp2);
    if(irc < 0) return irc;

    // add next edge as left-contour
    irc = MSH2D_brkl_bnd_add2 (brkls, if2, en2);
    if(irc < 0) return irc;

    irc = 0;
    goto L_exit;
  }

  // test if ip2o = right (-1) or left (1) or on (0) line ips-ipe
  spe = UT2D_sid_ptvc_tol (&pta[ip2o], brkls->ps, &brkls->vcl, &MSH2D_tol);
    // printf("  _sid_3pt-spo2 = %d p%d\n",spe,ip2o);


  //----------------------------------------------------------------
  // spe = 0 = ip2o is exact on line ips-ipe
  if(!spe) {
    // point ip2o is exact on edge. (eg ip2o is on straight line to ipe)
    ixTyp = 1; // point on ips-ipe, face right
    // set ipxf = point on ips-ipe
    ipxf = ip2o;
    // set if1,en1 = edge right of contour; endPt on line ips-ipe.
    if1 = if2;
    en1 = enp2;
    // set ixlf,ixle = edge right of contour; startPt on line ips-ipe.
    ixlf = if2;
    ixle = en2;
    // nothing to add to contours. Ignore point ip2o.
    // mark point iee unused
    MSH2D_pst_set (ip2o, 0);

    goto L_nxt;
  }


  //----------------------------------------------------------------
  // normal intersection.
  ixTyp = 0;
  if1 = if2;
 
  //..........................................
  if(spe < 0) {
    // ipo2 right of line; add enp2 to right contour; set en2 = ixf,ixe
    irc = MSH2D_brkl_bnd_add1 (brkls, if2, enp2);
    if(irc < 0) return irc;
    en1 = en2;

  //..........................................
  } else {
    // ipo2 left of line; add en2 to left contour; set enp2 = ixf,ixe
    irc = MSH2D_brkl_bnd_add2 (brkls, if2, en2);
    if(irc < 0) return irc;
    en1 = enp2;
  }

  goto L_nxt;


  //================================================================
  L_typ_1_:
  if(ixTyp != 1) goto L_typ_2_;
    printf(" L_typ_1_ -----------------------\n");
  // 1 = edge is parallel ips-ipe; face is right of line;
  // Using: if1,en1 = edge right of contour; start = cp1e, end = ipxf
  // Using: ipxf = point on ips-ipe (endPt of if1,en1)
  //   cp1e = lastPt-contour-right
  // add contour on right-side
    // printf(" if1=%d en1=%d ixlf=%d ixle=%d ipxf=%d\n",
            // if1,en1,ixlf,ixle,ipxf);

  // get if2 = neighbour-edge of if1,en1
  if2 = UFA_if_getNf_ife (if1, en1, MFNB__);

  // add face
  irc = MSH2D_brkl_bnd_addf (brkls, if2);
  if(irc < 0) return -1;

  // get en2 = edgeNr ending on ips-ipe
  en2 = UFA_esn_ipe (&MFAC__[if2], ipxf);
    // printf("  if2=%d en2=%d\n",if2,en2);

  // get prev & next edge of en2
  UFA_2esn_prvnxt (&enp2, &enn2, en2);

  // get ip2o = the opposite-point of if2,en2
  ip2o = UFA_ipo_fac_esn (&MFAC__[if2], enn2);

  // get spe = side of point ixps; 0=on, 1=left, -1=right
  spe = UT2D_sid_ptvc_tol (&pta[ip2o], brkls->ps, &brkls->vcl, &MSH2D_tol);
    // printf("  _sid_3pt-ip2s = %d p%d\n",spe,ip2o);

  //----------------------------------------------------
  // spe = 0 = exact on line ips-ipe
  if(!spe) {
    // 0 = ip2o and ipxf is on ips-ipe; edge parallel, but not end of line

    // check if ip2o = endPt of line
    if(ip2o == brkls->ipe) {
      // face if2 connects to ipe.
      // add prev edge to contour-right
      irc = MSH2D_brkl_bnd_add1 (brkls, if2, enp2);
      if(irc < 0) return irc;
      // set ipll = last point of contour-left
      ipll = ip2o;
      ixNxt = 0;  // finished
      goto L_typ_1l_;
    }

    // mark point ip2o unused
    MSH2D_pst_set (ip2o, 0);

    // get if3,en3 = nfb of if2,en2
    UFA_ife_getNf_ife (&if3, &en3, if2, en2, MFAC__, MFNB__);

    // add face (its on the left side !)
    irc = MSH2D_brkl_bnd_addf (brkls, if3);
    if(irc < 0) return -1;

    // get ipll = oppos. pt of if3,en3
    ipll = UFA_ipo_fac_esn (&MFAC__[if3], en3);
       // printf("  if3=%d en3=%d ipll=%d\n",if3,en3,ipll);

    ixTyp = 1; // pt on line
    ixNxt = 2; // fix left-bnd, continue with ixTyp=1;
    ipxf = ip2o;
    // set if1,en1 = edge right of contour; cp1e-ipxf
    if1 = if2;
    en1 = enp2;
    goto L_typ_1l_;

  //----------------------------------------------------
  } else if(spe < 0) {
    // -1 = ip2o is right of ips-ipe, ipxf is on ips-ipe.
    // add the prev. edge to right contour; keep ixTyp=1;
    // add prev edge to contour
    irc = MSH2D_brkl_bnd_add1 (brkls, if2, enp2);
    if(irc < 0) return irc;
    // set if1,en1 = last Edge from contourEndPt to pt on ips-ipe
    if1 = if2;
    en1 = en2;
    goto L_typ_1_; // redo check right side ...

  //----------------------------------------------------
  } else {
    // 1 = ip2o is left of ips-ipe, ipxf is on ips-ipe.
    // edge crosses ips-ipe. fix left contour,then continue with ixTyp=0
    ixTyp = 0;
    ixNxt = 1;           // continue with intersection
    // set ipll = last point of contour-left
    ipll = ip2o;
    // set ixf,ixe = intersecting-edge
    if1 = if2;
    en1 = enp2;
    goto L_typ_1l_;
  }


  //----------------------------------------------------------------
  // add contour on left-side until startpoint of contourEdge = ipll
  // Using: ixlf,ixle  last edge left of ips-ipe (cp2s-ipxf)
  // Using: ipll = find contour-left up to this point (if ixTyp=0)
  // Using: brkls->cp2s = startpoint of contour-left
  // Using: ipxf = point on ips-ipe
  // using: ixNxt=0: ipll=ipe; ixNxt=1: intersection; ..
  L_typ_1l_:
  // // get if1 = face with last edge of left-contour (has edge cp1e-ipxf)
  // if1 = brkls->cf2a[brkls->c2Nr - 1];
  // // en1 = brkls->ce2a[brkls->c2Nr - 1];
    // printf("L_typ_1l: ----------------------------------\n");
    // printf(" ixTyp=%d ixNxt=%d ixlf=%d ixle=%d ipxf=%d cp2s=%d ipll=%d\n",
             // ixTyp,   ixNxt,   ixlf,   ixle,   ipxf,brkls->cp2s,ipll);
    // if((ipll==195)&&(brkls->c1Nr>10)) return -99;


  if(brkls->cp2s == ipll) {
    // end reached.

    // test if ipxf reached.
    // get ip2s = startPt of ixlf,ixle
    ip2s = UFA_ips_fac_esn (&MFAC__[ixlf], ixle);
    if(ip2s != ipxf) {
      // get nfb and its prev;edge. set this edge as ixlf,e
      irc = UFA_ife_getNf_ife (&if2, &enp2, ixlf, ixle,  MFAC__, MFNB__);
        // printf(" ip2s=%d if2=%d &enp2=%d irc=%d\n",ip2s,if2,enp2,irc);
      ixlf = if2;
      ixle = UFA_esn_prv (enp2);
    }

    // exit if (ipll == ipe)
    if(ixNxt == 0) goto L_exit;

    // ixNxt=1: next face has intersection, continue with new ixTyp=0
    // if1,en1 = intersecting-edge
    if(ixNxt == 1) goto L_nxt;

    // ixNxt=2: next point is also on ips-ipe;
    // keep ipxf = point on ips-ipe
    // keep ixlf,ixle = edge left of contour; ipxf-cp2s
    if(ixNxt == 2) goto L_typ_1_;

TX_Error("MSH2D_brkl_bnd_2_ ixTyp=%d-L_typ_1l5:E1",ixTyp); return -99;
  }


  // get if2 = neighbour-face of ixlf,ixle
  // if2 = UFA_if_getNf_ife (ixlf, ixle, MFNB__);
  irc = UFA_ife_getNf_ife (&if2, &enp2, ixlf, ixle,  MFAC__, MFNB__);
  if(irc < 0) {TX_Error("MSH2D_brkl_bnd_2_ E5"); return -99;}

  // get prev & nxt
  UFA_2esn_prvnxt (&enn2, &en2, enp2);
    // printf(" if2=%d enp2=%d en2=%d enn2=%d\n",if2,enp2,en2,enn2);

  // add face
  irc = MSH2D_brkl_bnd_addf (brkls, if2);
  if(irc < 0) return -1;

  // add next edge to contour;
  irc = MSH2D_brkl_bnd_add2 (brkls, if2, enn2);
  if(irc < 0) return irc;

  ixlf = if2;
  ixle = en2;
  goto L_typ_1l_;


  //================================================================
  L_typ_2_:
    // 2 = edge is parallel ips-ipe; face is left of line;
    // endpoint of edge ixf,ief is on line ips-ipe.
    // find next face on left side
      // printf(" L_typ_2_: ixf=%d ixe=%d cp1e=%d\n",
              // brkls->ixf,brkls->ixe,brkls->cp1e);
    if1 = brkls->ixf;
    en1 = UFA_esn_nxt (brkls->ixe);
      // printf(" L_typ_2_: if1=%d en1=%d\n",if1,en1);
    // get neighbour-face
    if2 = UFA_if_getNf_ife (if1, en1, MFNB__);
    // get en2 = edgeNr starting with cp2s
    en2 = UFA_esn_ips (&MFAC__[if2], brkls->cp2s);
      // printf(" L_typ_2_: if2=%d en2=%d\n",if2,en2);
    // get ip2o = opposite point
    UFA_ipOpp_fac_esn (&ip2o, &i1, &MFAC__[if2], en2);
      // printf(" ip2o=%d\n",ip2o);
    // test if ipo2 == ipe
    if(ip2o == brkls->ipe) {
      // add face
      irc = MSH2D_brkl_bnd_addf (brkls, if2);
      if(irc < 0) return -1;
      // add prev to left -contour
      enp2 = UFA_esn_prv (en2);
      irc = MSH2D_brkl_bnd_add2 (brkls, if2, enp2);
      if(irc < 0) return irc;
      ixNxt = 0;
      goto L_exit;
    }

    // test if ip2o = right (-1) or left (1) or on (0) line ips-ipe
    spe = UT2D_sid_ptvc_tol (&pta[ip2o], brkls->ps, &brkls->vcl, &MSH2D_tol);
      // printf("  _sid_3pt-spo2 = %d p%d\n",spe,ip2o);
    if(spe == 0) {
      // edge starts also on ips-ipe; add to faceList, ignore;
      // add face
      irc = MSH2D_brkl_bnd_addf (brkls, if2);
      if(irc < 0) return -1;

// TODO: point unused - ip2o can be deleted ?
    TX_Print("MSH2D_brkl_bnd_2_ I3-%d",ip2o);

      brkls->ixf = if2;
      brkls->ixe = UFA_esn_nxt (en2);
      brkls->cp1e = ip2o;
      goto L_typ_2_;

    } else if(spe == 1) {
      // edge has only one point at ips-ipe; add to faceList & to left-contour;
      // add face
      irc = MSH2D_brkl_bnd_addf (brkls, if2);
      if(irc < 0) return -1;
      // add prev edge to left-boundary
      enp2 = UFA_esn_prv (en2);
      irc = MSH2D_brkl_bnd_add2 (brkls, if2, enp2);
      if(irc < 0) return irc;

      brkls->ixf = if2;
      brkls->ixe = en2;
      brkls->cp2s = ip2o;
      goto L_typ_2_;

    } else {
      TX_Error("MSH2D_brkl_bnd_2_ ixTyp=%d",ixTyp); return -99;
    }


  //================================================================
  L_exit:
    // TESTBLOCK
    // printf("ex-MSH2D_brkl_bnd_2_ irc=%d ips=%d\n",irc,brkls->ips);
    // if(brkls->ips == 0) {
      // MSH2D_dump_brkls (brkls, "sb");
      // return MSG_STD_ERR (ERR_internal, "exit; irc=%d ips=%d",irc,brkls->ips);
    // }
    // END TESTBLOCK

  return irc;

}


//================================================================
  int MSH2D_brkl_bnd_23 (brklStru *brkls) {
//================================================================
// faces only right of ips-ipe
// contours done up to ipxf on right side and cp2s on left side.
// follow ips-ipe until ipe.
// startPt ipxf is on line ips-ipe;
// Intput:
//   ipxf        startPoint of ixlf,ixle; 
//   ixlf,ixle   edge ipxf-cp2s
// Output: 
//   RetCod:
//      0       OK; contours & faces in brkls 

 
  int     irc, if2, en2, enn2, enp2, ip2o, spe;


  // printf("MSH2D_brkl_bnd_23 \n");


  L_nxt:
  // get nfb
  irc = UFA_ife_getNf_ife (&if2, &enn2, brkls->ixf, brkls->ixe, MFAC__, MFNB__);
  if(irc < 0) {TX_Error("MSH2D_brkl_bnd_23 E1"); return -99;}
    // printf(" if2=%d enp2=%d\n",if2,enp2);

  // add face
  irc = MSH2D_brkl_bnd_addf (brkls, if2);
  if(irc < 0) return -1;

  // get ipo2 = opposite point and en2 = prev enn2
  UFA_ipOpp_fac_esn (&ip2o, &en2, &MFAC__[if2], enn2);
    // printf(" ip2o=%d en2=%d\n",ip2o,en2);

  // get en2 = prev enn2
  enp2 = UFA_esn_prv (en2);

  // test endpoint reached. Yes: add to contour-left & -right;
  if(ip2o == brkls->ipe) {
    // add edge & next to contours
    irc = MSH2D_brkl_bnd_add1 (brkls, if2, enp2);
    if(irc < 0) return irc;
    // irc = MSH2D_brkl_bnd_add2 (brkls, if2, en2);
    // if(irc < 0) return irc;
    brkls->ixTyp = 16;
    irc = 0;
    goto L_exit;
  }

  // get spe = side of point ixps; 0=on, 1=left, -1=right
  spe = UT2D_sid_ptvc_tol (&brkls->pta[ip2o], brkls->ps, &brkls->vcl, &MSH2D_tol);
    // printf("  _sid_3pt-ip2s = %d p%d\n",spe,ip2o);

  // -1 = ipo2 = right of ips-ipe: add enp2 to contour-right; continue L_nxt;
  if(spe < 0) {
    irc = MSH2D_brkl_bnd_add1 (brkls, if2, enp2);
    if(irc < 0) return irc;
    brkls->ixf =  if2;
    brkls->ixe =  en2;
    goto L_nxt;
  }

  // // 0 = ipo2 on line: add to contour-left; set ixf,ixe; continue L_nxt;
  // irc = MSH2D_brkl_bnd_add2 (brkls, if2, en2);
  if(irc < 0) return irc;

  brkls->ixf =  if2;
  brkls->ixe =  enp2;
  // brkls->ipxf = ip2o;
  goto L_nxt;



  //----------------------------------------------------------------
  L_exit:
    // TESTBLOCK
    // printf("ex-MSH2D_brkl_bnd_23 irc=%d\n",irc);
    // MSH2D_dump_brkls (brkls, "ex-MSH2D_brkl_bnd_23");
    // if(brkls->ips == 0)
    // return -99;
    // END TESTBLOCK

  return 0;

}


//================================================================
  int MSH2D_brkl_bnd_24 (brklStru *brkls) {
//================================================================
// faces only left of ips-ipe
// contours done up to cp1e on right side and ipxf on left side.
// follow ips-ipe until ipe.
// startPt ipxf is on line ips-ipe;
// Intput:
//   ipxf        startPoint of ixf,ixe; 
//   ixf,ixe     edge ipxf-cp1e
// Output:
//   RetCod:
//      0       OK; contours & faces in brkls



  int     irc, if2, en2, enn2, enp2, ip2o, spe;


  // printf("MSH2D_brkl_bnd_24 \n");


  L_nxt:
  // get nfb
  irc = UFA_ife_getNf_ife (&if2,&enp2, brkls->ixlf,brkls->ixle, MFAC__, MFNB__);
  if(irc < 0) {TX_Error("MSH2D_brkl_bnd_24 E1"); return -99;}
    // printf(" if2=%d enp2=%d\n",if2,enp2);

  // add face
  irc = MSH2D_brkl_bnd_addf (brkls, if2);
  if(irc < 0) return -1;

  // get ipo2 = opposite point and enn2 = prev enp2
  UFA_ipOpp_fac_esn (&ip2o, &enn2, &MFAC__[if2], enp2);
    // printf(" ip2o=%d\n",ip2o);

  // get en2 = prev enn2
  en2 = UFA_esn_prv (enn2);

  // test endpoint reached. Yes: add to contour-left & -right;
  if(ip2o == brkls->ipe) {
    // add edge & next to contours
    // irc = MSH2D_brkl_bnd_add1 (brkls, if2, en2);
    // if(irc < 0) return irc;
    irc = MSH2D_brkl_bnd_add2 (brkls, if2, enn2);
    if(irc < 0) return irc;
    brkls->ixTyp = 16;
    irc = 0;
    goto L_exit;
  }

  // get spe = side of point ixps; 0=on, 1=left, -1=right
  spe = UT2D_sid_ptvc_tol (&brkls->pta[ip2o], brkls->ps, &brkls->vcl, &MSH2D_tol);
    // printf("  _sid_3pt-ip2s = %d p%d\n",spe,ip2o);

  // 1 = ipo2 = left of ips-ipe: add enn2 to contour-left; continue L_nxt;
  if(spe > 0) {
    irc = MSH2D_brkl_bnd_add2 (brkls, if2, enn2);
    if(irc < 0) return irc;
    brkls->ixlf =  if2;
    brkls->ixle =  en2;
    goto L_nxt;
  }

  // // 0 = ipo2 on line: add to contour-right; set ixlf,ixle; continue L_nxt;
  // irc = MSH2D_brkl_bnd_add1 (brkls, if2, en2);
  // if(irc < 0) return irc;

  brkls->ixlf =  if2;
  brkls->ixle =  enn2;
  // ipxf = ip2o;
  goto L_nxt;



  //----------------------------------------------------------------
  L_exit:
    // TESTBLOCK
    // printf("ex-MSH2D_brkl_bnd_24 irc=%d\n",irc);
    // MSH2D_dump_brkls (brkls, "ex-MSH2D_brkl_bnd_24"); 
    // if(brkls->ips == 0)
    // return -99;
    // END TESTBLOCK

  return 0;

}


//======================================================================
  int MSH2D_brkl_bnd_add1 (brklStru *brkls, int if1, int en1) {
//======================================================================
// add edge to right-contour; set brkls->cp1e

  int    ii;

  ii = brkls->c1Nr;

  // printf("MSH2D_brkl_bnd_add1 if1=%d en1=%d c1Nr=%d\n",if1,en1,ii);

  if(if1 < 0) {TX_Error("MSH2D_brkl_bnd_add1 E1"); return -91;}

  brkls->cf1a[ii] = if1;
  brkls->ce1a[ii] = en1;

  brkls->c1Nr += 1;
	if(brkls->c1Nr >= brkls->siz) {
    TX_Error("MSH2D_brkl_bnd_add1 E1");
    return -1;
  }

  // set brkls->cp1e = endPt of right-contour
  brkls->cp1e = UFA_ipe_fac_esn (&MFAC__[if1], en1);

  return 0;

}

//======================================================================
  int MSH2D_brkl_bnd_add2 (brklStru *brkls, int if1, int en1) {
//======================================================================
// add edge to left-contour; set brkls->cp2s

  int    ii;

  ii = brkls->c2Nr;

  // printf("MSH2D_brkl_bnd_add2 if1=%d en1=%d c1Nr=%d\n",if1,en1,ii);
  // if(if1 == 344) return -99;

  if(if1 < 0) {TX_Error("MSH2D_brkl_bnd_add2 E1"); return -91;}

  brkls->cf2a[ii] = if1;
  brkls->ce2a[ii] = en1;

  brkls->c2Nr += 1;
  if(brkls->c1Nr >= brkls->siz) {
    TX_Error("MSH2D_brkl_bnd_add2 E1");
    return -1;
  }

  // set brkls->cp2s = startPt of left-contour
  brkls->cp2s = UFA_ips_fac_esn (&MFAC__[if1], en1);

  return 0;

}


//======================================================================
  int MSH2D_brkl_bnd_addf (brklStru *brkls, int if1) {
//======================================================================
// add to facelist (uniq)
// if eCxa=NULL then init

  int    irc, i1, ii, *ia;

  // printf("MSH2D_brkl_bnd_addf if1=%d\n",if1);


  ii = brkls->fiNr;
  ia = brkls->fia;

  for(i1=0; i1 < ii; ++i1) if(ia[i1] == if1) return 1;

  ia[ii] = if1;

  brkls->fiNr += 1;

  if(brkls->fiNr >= brkls->siz) {
    TX_Error("MSH2D_brkl_bnd_addf E1");
    return -1;
  }

  return 0;

}



//======================================================================
  int MSH2D_brkl_flip (int if1, int en1, int ips, int ipe, int edgTyp) {
//======================================================================


  int      if2, en2;

  // get neigbour-face and common edge of if1,en1
  UFA_ife_getNf_ife (&if2, &en2, if1, en1, MFAC__, MFNB__);
    // printf(" if2=%d en2=%d\n",if2,en2);


  // flip
  UFA_2fac_flip (if1, en1, if2, en2, MFAC__, MFNB__);

  en1 = UFA_esn_ips (&MFAC__[if1], ips); // get edge# of ips in ifc1
  en2 = UFA_esn_ips (&MFAC__[if2], ipe); // get edge# of ipe in ifc2

  UFA_fnb_set_if_es (MFNB__, if1, en1, edgTyp);
  UFA_fnb_set_if_es (MFNB__, if2, en2, edgTyp);


  return 0;

}
 

//=====================================================================
  int MSH2D_brkl_set (int if1, int en1, int ips, int ipe, int edgTyp) {
//=====================================================================
// set edgType of edge (fnb of the common line)
 
  int      if2, en2;


  // get neigbour-edge
  UFA_ife_getNf_ife (&if2, &en2, if1, en1, MFAC__, MFNB__);
    // printf(" _brkl_set-if2=%d en2=%d\n",if2,en2);

  UFA_fnb_set_if_es (MFNB__, if1, en1, edgTyp);


  if(if2 >= 0)
    UFA_fnb_set_if_es (MFNB__, if2, en2, edgTyp);


  return 0;

}


//===========================================================================
  int MSH2D_brkl_cpy (int *fan, MemTab(Fac3) *mfi) {
//===========================================================================
// add mfi = new faces to mfo (overwriting deleted faces)
// Output:
//   fan      list of new faces

// TODO: deleted fac+fnb records in eine Liste eintragen;
// bei add-new-records diese Liste zuerst benutzen; wenn leer normales add.



  int     iii, iie, i1, i2, fni, fNr;
  Fac3    *fai;


  fai = MEMTAB_DAT(mfi);
  fni = MEMTAB_IND(mfi);

  fNr = MEMTAB_IND (&MSH2D_ts->fac);

  // printf("MSH2D_brkl_cpy fni=%d fNr=%d\n", fni, fNr);

  iii = 0;

  //----------------------------------------------------------------
  // overwrite deleted faces (with .st = "quadIndex-unknown")
  // find next deleted face
  // fao = MEMTAB_DAT(mfo);
  for(i1=0; i1<fNr; ++i1) {
    if(MFAC__[i1].st >= 0) continue;   // skip active faces

    // MFAC__[i1] = deleted face; overwrite with fai[iii]
    MFAC__[i1] = fai[iii];
    MFAC__[i1].st = MSH2D_qNr; // "quadIndex-unknown"
      // printf(" brkl_cpy %d -> %d\n",iii,i1);
    // invalidate fnb
    // MFNB__[i1] = UFA_NUL;
    MFNB__[i1].i1 = MSH_EDGLN_ERR;
    MFNB__[i1].i2 = MSH_EDGLN_ERR;
    MFNB__[i1].i3 = MSH_EDGLN_ERR;
    // add ifac to list fan
    fan[iii] = i1;
    // next
    ++iii;
    if(iii >= fni) break;
  }
    // printf(" iii=%d fni=%d\n",iii,fni);

  //----------------------------------------------------------------
  if(iii >= fni) {
    // all new faces have been copied.
    // remove all deleted faces at end of faceTab
      // printf(" remove i1=%d fNr=%d\n",i1,fNr);
    for(i2 = (fNr - 1); i2 >= iii; --i2) {
      if(MFAC__[i2].st >= 0) break;   // skip active faces
        // printf(" remove i2=%d\n",i2);
      MSH2D_ts->fac.rNr = i2; // remove faces at end of faceTab
    }
    goto L_exit;
  }



  //----------------------------------------------------------------
  // add remaining faces (with .st = "quadIndex-unknown")
  if(iii <= fNr) {
    // fni = nr of records to store;
    // iii = nr of records already stored in MFAC__;

    // check space for additional fac's, fnb's (fnb must have same size as fac)
    // get free space in MFAC__
    if(MEMTAB_RFREE(&MSH2D_ts->fac) < iii) {
      // MSH2D_spc_ck1 (iii);
      TX_Error("MSH2D_brkl_cpy E1");
      return -99;
    }


    for(i1=iii; i1<fni; ++i1) {
      // add face
      UFA_add_fac_st (&MSH2D_ts->fac,
                      fai[i1].i1, fai[i1].i2, fai[i1].i3, MSH2D_qNr);
      // add ifac to list fan
      fan[i1] = MSH2D_ts->fac.rNr - 1;  // MEMTAB_IND(mfo) - 1;
        // printf(" add-mfo fan[%d] = %d\n",i1,fan[i1]);
    }
  }



  //================================================================
  L_exit:


    // TESTBLOCK
    // MSH2D_dump__ ("ex-MSH2D_brkl_cpy");  MAKES STOP !
    // END TESTBLOCK


  return 0;

}


//================================================================
  int MSH2D_quad_set_mod (int iqx) {
//================================================================
// MSH2D_quad_set_mod         set quad and its faces to modified
// see also MSH2D_quad_nf_mod MSH2D_quad_set_mf

  int   i1,  ii, fNr;

  // printf("MSH2D_quad_set_mod %d\n",iqx);

  ii = 0;
  fNr = MEMTAB_IND (&MSH2D_ts->fac);

  for(i1=0; i1<fNr; ++i1) {
    if(MFAC__[i1].st != iqx) continue; // skip value "quadIndex-unknown"
    // found face in quad iqx; set to modified
    MFAC__[i1].st = MSH2D_qNr;
    ++ii;
      // printf(" quad_set_mod-f%d mod; %d\n",i1,ii);
    if(ii >= 2) break;
  }

  MSH2D_qma[iqx] = 0;  // set modified

  return 0;

}

 
//================================================================
  int MSH2D_quad_nf_mod (int cfNr, int *cfa) {
//================================================================
// set all quads of faces to delete to modified !

  int      i1, iqx;

  // printf("MSH2D_quad_nf_mod %d\n",cfNr);
// if(cfNr==5) return -99;

  for(i1=0; i1<cfNr; ++i1) {
    iqx = MFAC__[cfa[i1]].st;
    if(iqx == MSH2D_qNr) continue; // skip value "quadIndex-unknown"
      // printf(" qma_mod %d cfa=%d iqx=%d qma=%d\n",i1,cfa[i1],iqx,MSH2D_qma[iqx]);
    // quad and its faces to modified
    MSH2D_quad_set_mod (iqx);
  }


  return 0;

}


//================================================================
  int MSH2D_del__ () {
//================================================================
// MSH2D_del__    delete unused gridpoints
// remove all unused points of pa2; update .fac;
// DOES NOT FIX .tab and .ipa;
 

  int      irc, i1, itNr, ib, ptNr, *ptgStat, fNr, ipx;
  int      ipSta, ipTo, ipFrom, ipEnd, ipNr, *ipa, ipLast;
  IndTab   *itTab;
  Point2   *pta;
  // Fac3     *fa;


  // printf("MSH2D_del__ \n");
  // UFA_fnb_dump__ (MFAC__, MFNB__,
                  // MEMTAB_IND(&MSH2D_ts->fac),  "_del__");


  itNr  = MEMTAB_IND (&MSH2D_ts->tab);  // nr of contours
  itTab = MEMTAB_DAT (&MSH2D_ts->tab);  // contours

  fNr = MEMTAB_IND (&MSH2D_ts->fac);  // ind of first face
  // fa  = MEMTAB_DAT (&MSH2D_ts->fac);

  ptNr = MEMTAB_IND (&MSH2D_ts->pa2);
  pta  = MEMTAB_DAT (&MSH2D_ts->pa2);

  ipa = MEMTAB_DAT (&MSH2D_ts->ipa);


  // find MSH_GRIDBOX-record
    for(ib=0; ib<itNr; ++ib) {
    if(itTab[ib].typi != MSH_GRIDBOX) continue;

    ipSta = itTab[ib].ibeg;
    ipNr  = itTab[ib].iNr;
    ipEnd = ipSta + ipNr - 1;
    goto L_1;
  }
  return 1;    // no MSH_GRIDBOX record in .tab


  //----------------------------------------------------------------
  // set ipSta = index first gridpoint
  L_1:
	  // printf(" ipSta=%d ipEnd=%d ipNr=%d ptNr=%d\n",ipSta,ipEnd,ipNr,ptNr);


  // get ptgStat = temp-space; one int for one gridpoint
  ptgStat = (int*) MEM_alloc_tmp (ipNr * (int)sizeof(int));

  // init ptgStat to -2 (unused)
  for(i1=0; i1<ipNr; ++i1) ptgStat[i1] = -2;


  // loop tru all faces; set ptgStat for used points to -1 (used);
  for(i1=0; i1<fNr; ++i1) {
    ipx = MFAC__[i1].i1 - ipSta; if(ipx >= 0) ptgStat[ipx] = -1;
    ipx = MFAC__[i1].i2 - ipSta; if(ipx >= 0) ptgStat[ipx] = -1;
    ipx = MFAC__[i1].i3 - ipSta; if(ipx >= 0) ptgStat[ipx] = -1;
  }
    // TESTBLOCK
    // for(i1=0; i1<ipNr; ++i1)
       // printf(" ptgStat[%d]=%d / %d\n",i1,ptgStat[i1],i1+ipSta);
    // END TESTBLOCK


  // loop tru ptgStat; delete unused point (copy last used over unused)
  // -2 unused; -1 used; >=0 new-index;
  ipTo = -1;          // ipSta - 1;
  ipFrom = ipNr - 1;  // ipEnd;
  goto L_59;


 
  //----------------------------------------------------------------
  L_50:
    // find ipTo = next-unused-point (ptgStat=0);
      // printf(" ipTo=%d\n",ipTo);


    // find ipFrom = last-used-point (-1)
  L_51:
    if(ptgStat[ipFrom] != -1) {
      --ipFrom;
      if(ipFrom <= ipTo) goto L_up_0;
      goto L_51;
    }
      // printf(" ipFrom=%d\n",ipFrom);

    // copy point; overwrite next-unused-point with last-used-point
    pta[ipTo + ipSta] =  pta[ipFrom + ipSta];
      // printf(" copy %d = %d\n",ipTo+ipSta,ipFrom+ipSta);

    // store new pointindex in ptgStat
    ptgStat[ipFrom] = ipTo;

    // set ipTo=used (-1)
    ptgStat[ipTo] = -1;

    // update ipa
    ipa[ipSta + ipTo] = ipSta + ipTo;
 

    // find ipTo = next-unused-point (ptgStat=0);
  L_59:
    ++ipTo;
    if(ipTo >= ipNr) goto L_up_0;
    if(ptgStat[ipTo] != -2) goto L_59;  // -2 = unused = ignore

    goto L_50;


  //----------------------------------------------------------------
  // loop tru all faces; update points
  //   if ptgStat[i1|i2|i3] >= 0) then i1|i2|i3 = ptgStat[i1|i2|i3]
  L_up_0:
    for(i1=0; i1<fNr; ++i1) {
      ipx = MFAC__[i1].i1 - ipSta;
      if(ipx >= 0) {
        if(ptgStat[ipx] >= 0)    MFAC__[i1].i1 = ptgStat[ipx] + ipSta;
      }

      ipx = MFAC__[i1].i2 - ipSta;
      if(ipx >= 0) {
        if(ptgStat[ipx] >= 0)    MFAC__[i1].i2 = ptgStat[ipx] + ipSta;
      }

      ipx = MFAC__[i1].i3 - ipSta;
      if(ipx >= 0) {
        if(ptgStat[ipx] >= 0)    MFAC__[i1].i3 = ptgStat[ipx] + ipSta;
      }

    }

    // find ipLast = last used point
    ipLast = -1;
    for(i1=ipNr-1; i1>=0; --i1) {
      if(ptgStat[i1] != -1) continue;
      ipLast = i1;
      break;
    }
      // printf(" ipLast = %d\n",ipLast);


    // update MSH_GRIDBOX-record
    itTab[ib].iNr = ipLast + 1;

    // update nr of 2D-point. ipTo is index of last-used-point;
    ipLast = ipLast + ipSta + 1;
    MSH2D_ts->pa2.rNr = ipLast;

    // remove unused ipas
    MSH2D_ts->ipa.rNr = ipLast;



  //----------------------------------------------------------------
    // TESTBLOCK   ex MSH2D_del__
    // printf("L_up_0: \n");
    // for(i1=0; i1<ipNr; ++i1) printf(" ptgStat[%d]=%d\n",i1,ptgStat[i1]);
    // UFA_fnb_dump__ (MFAC__, MFNB__,
                    // MEMTAB_IND(&MSH2D_ts->fac),  "ex-_del__");
    // MemTab_dump (&MSH2D_ts->ipa, "ex-_del__");
    // END TESTBLOCK


  return 0;

}


//================================================================
  int MSH2D_fdel_bnd (int cNr, int *cfa, char *cea) {
//================================================================
// delete all faces inside closed contour (set face.st = -1)
// Input:
//   cNr    nr of edges cfa,cea
// set fac.st = -1
// 

  int      irc, i1, i2, i3, ii, f1, e1, fn, en;
  TessStru *ts1 = MSH2D_ts;
  // Fac3     *fTab;


  // printf("MSH2D_fdel_bnd %d\n",cNr);
  // for(i1=0;i1<cNr;++i1)printf(" edge-%d f-%d e-%d\n",i1,cfa[i1],cea[i1]);


  // fTab = MEMTAB_DAT (&ts1->fac);

  // delete all faces of the edges of contour
  // loop tru contour
  for(i1=0; i1<cNr; ++i1) {
    f1 = cfa[i1];
    e1 = cea[i1];
      // printf(" nxt-cSeg %d %d\n",f1,e1);


    // delete this face f1
    L_f_del:
    MFAC__[f1].st = -1;

    // loop tru next 2 edges of face
    L_edg__:
    for(i2=0; i2<2; ++i2) {
      // set e1 = next edge
      e1 = UFA_esn_nxt(e1);
      // check if edge is also part of the boundary;
      for(i3=0; i3<cNr; ++i3) {
        if(cfa[i3] != f1) continue;   
        if(cea[i3] != e1) continue;   
        // f1,e1 is part of the contour; stop here (do not delete)
        goto L_edg_nxt;
      }
      // Get fnb fn,en of this edge
      irc = UFA_ife_getNf_ife (&fn, &en, f1,e1, MFAC__, MFNB__);
      if(irc){printf("MSH2D_fdel_bnd I1");return -1;}
      // if(irc){TX_Error("MSH2D_fdel_bnd E1");return -1;}
        // printf("getNf_ife-irc=%d fn=%d en=%d st=%d\n",irc,fn,en,MFAC__[fn].st);

      // if face is already deleted: continue.
      if(MFAC__[fn].st < 0) continue;
      f1 = fn;
      e1 = en;
      // delete face
      goto L_f_del;

      L_edg_nxt:
      continue;
    }

    L_nxt_fac:
    continue;
  }


  return 0;

}

 
//======================================================================
  int MSH2D_ctess_e (MemTab(Fac3) *fca, int cfNr, int *cfa, char *ena) {
//======================================================================
// get contour = points of edges ips-ipe and tesselate contour
// Input:
//   cfNr       nr of edges
//   cfa,ena    contour of edges, closing edge not included.


  int      irc, i1, i2, ipNr, *ipa;
  Point2   *pta;
  Fac3     *fa;
  TessStru *ts1 = MSH2D_ts;


  printf("MSH2D_ctess_e %d\n",cfNr);


  fa = MEMTAB_DAT (&ts1->fac);
  pta = MEMTAB_DAT (&ts1->pa2);
  

  // get space for points
  ipNr = cfNr + 1;
  ipa = (int*) MEM_alloc_tmp (ipNr * (int)sizeof(int));
  if(ipa == NULL) {TX_Error("MSH2D_ctess_e EOM-1"); return -1;}

  // add startpoint of all edges
  for(i1=0; i1<cfNr; ++i1) {
    ipa[i1] = UFA_ips_fac_esn (&fa[cfa[i1]], ena[i1]);
  }

  // add endpoint of last edge
  i1 = cfNr - 1;
  ipa[cfNr] = UFA_ipe_fac_esn (&fa[cfa[i1]], ena[i1]);

    // for(i1=0;i1<ipNr;++i1) printf(" ipa[%d] = %d\n",i1,ipa[i1]);

  // tesselate
  irc = MSH2D_ctess__ (fca, ipa, ipNr, pta);

 // set .st of all faces to modified = MSH2D_qNr
  i2 = MEMTAB_IND(fca);
  for(i1=0; i1<i2; ++i1) fca->data[i1].st = MSH2D_qNr;


    // for(i1=0; i1<i2; ++i1) UFA_fac_dump1 (&fca->data[i1], "fca%d",i1);


  return 0;

}


//================================================================
  int MSH2D_view_inpt (int ptNr, int *ipa, Point2 *pa2) {
//================================================================
// MSH2D_view_cv_inpt             display curve of n indexed 2D-points

  int    i1;
  double zVal = 0.01;
  // Point2 *pa2;
  Point  *pta;


  printf("MSH2D_view_inpt %d\n",ptNr);

  // pa2 = MEMTAB_DAT (&MSH2D_ts->pa2);

  // get space for inptNr 3D-points
  pta = (Point*) MEM_alloc_tmp (ptNr * (int)sizeof(Point));

  // UT3D_npt_inpt2
  for(i1=0; i1<ptNr; ++i1) {
    // memcpy (&pta[i1], &pt2a[inpta[i1]], sizeof(Point2));
    // pta[i1].z = zVal;
    UT3D_pt_ptz (&pta[i1], &pa2[ipa[i1]], zVal);
  }


  GL_view_npt (pta, ptNr, ATT_PT_GREEN, -1L);
  // GL_view_nsymB  (ptax, ptNr, ATT_COL_RED, SYM_STAR_S, -1L);


  return 0;

}


//================================================================
  int MSH2D_view_cv_inpt (int inptNr, int *inpta) {
//================================================================
// MSH2D_view_cv_inpt             display curve of n indexed 2D-points

  int    i1;
  double zVal = 0.;
  Point2 *pt2a;
  Point  *pta;


  printf("MSH2D_view_cv_inpt %d\n",inptNr);

  pt2a = MEMTAB_DAT (&MSH2D_ts->pa2);

  pta = (Point*) MEM_alloc_tmp (inptNr * (int)sizeof(Point));

  // UT3D_npt_inpt2
  for(i1=0; i1<inptNr; ++i1) {
    memcpy (&pta[i1], &pt2a[inpta[i1]], sizeof(Point2));
    pta[i1].z = zVal;
  }


  GR_Disp_cv (pta, inptNr, 7); // Typ_Att_top1);
  // GL_DrawPoly2D

  // GR_Disp_pTab (inptNr, pta, SYM_TRI_B, 11); //ATT_COL_BLUE);

  return 0;

}
 

//===============================================================================
  int MSH2D_view_nfac (TessStru *ts1, char *opts, int if1, int fNr, char *info) {
//===============================================================================
// display faces, points (numbers) ..
//   if1      index of first face to display
//   fNr      nr of faces to display or -1 = all
//   opts v   disp normalvector in gravity-centerPoint
//        f   disp faceNr in gravity-centerPoint
//        b   disp bounds of faces
//        s   disp triangle shaded
//        P   disp points as (SYM_STAR_S,ATT_COL_BLUE)
//        n   disp point-numbers

  int      i1, pNr;
  Fac3     *fa;
  Point2   *pa2;
  Point    *pa3;		


  if(info) printf("MSH2D_view_nfac |%s| %s\n",opts,info);

  if(!ts1) return -1;

  // MSH2D_view_nfac (&ts1->fac, if1, fNr, &ts1->pa2, opts);

  fa = MEMTAB_DAT (&ts1->fac);
  if(fNr < 0) fNr = MEMTAB_IND (&ts1->fac);

  pa2 = MEMTAB_DAT (&ts1->pa2);
  pNr = MEMTAB_IND (&ts1->pa2);


  // make 3D-points from 2D-points
  pa3 = (Point*) MEM_alloc_tmp (pNr * (int)sizeof(Point));
  for(i1=0; i1<pNr; ++i1) pa3[i1] = UT3D_pt_pt2 (&pa2[i1]);


  UFA_view__ (fa, if1, fNr, pa3, pNr, opts, -1L);


  return 0;

}

 
//================================================================
  int MSH2D_dump__ (char *inf) {
//================================================================
// MSH2D_dump__              UFA_fnb_dump__ & UFA_nfb_ck__

  int   i1, i2, fNr;


  fNr = MEMTAB_IND(&MSH2D_ts->fac);

  // dump faces and neighbours
  UFA_fnb_dump__ (MFAC__, MFNB__, fNr, inf);

  // check faces for zero-size (linear face)
  i1 = UFA_fac_ck_zero_2D (MFAC__, fNr,
                           MEMTAB_DAT (&MSH2D_ts->pa2), &MSH2D_tol, inf);

  // check neighbours
  i2 = UFA_nfb_ck__ (MFAC__, MFNB__, fNr, inf);
  if((i1)||(i2)) return -99;

  return 0;

}


//================================================================
  int MSH2D_dump_qma () {
//================================================================

  int    i1, *ia;

  printf("MSH2D_dump_qma %d\n",MSH2D_qNr);
  for(i1=0; i1<MSH2D_qNr; ++i1) printf(" qma[%d] = %d\n",i1,MSH2D_qma[i1]);

  return 0;

}
 

//================================================================
  int MSH2D_dump_brkls (brklStru *brkls, char *opts) {
//================================================================
// MSH2D_dump_brkls       dump crossing-edges
// opts: s display brkls.fia as faces
//       b display brkls.cf1a & .cf2a as curves
// TODO: use UFA_Disp_nEdg2 for display curves ?

  int      i1, i2, *ia, igl;
  double   zVal=0.01;
  Point2   pta[2];

  printf("MSH2D_dump_brkls ips=%d opts=|%s|\n",brkls->ips,opts);



  //----------------------------------------------------------------
  // mode 0
  for(i1=0; i1 < brkls->c1Nr; ++i1) {
    printf(" cf1a[%d] ifc=%d en=%d\n",i1,
            brkls->cf1a[i1], brkls->ce1a[i1]);
  }

  for(i1=0; i1 < brkls->c2Nr; ++i1) {
    printf(" cf2a[%d] ifc=%d en=%d\n",i1,
            brkls->cf2a[i1], brkls->ce2a[i1]);
  }

  for(i1=0; i1 < brkls->fiNr; ++i1) {
    printf(" fia[%d] ifc=%d\n",i1,
            brkls->fia[i1]);
  }

  printf(" ixTyp=%d ixf=%d ixe=%d\n",brkls->ixTyp,brkls->ixf,brkls->ixe);
  printf(" cp1e=%d cp2s=%d ipxf=%d\n",brkls->cp1e,brkls->cp2s,brkls->ipxf);
  printf("   ixlf=%d ixle=%d\n",brkls->ixlf,brkls->ixle);


  //----------------------------------------------------------------
  if(strlen(opts) < 1) goto L_end;

  GL_view_ini__ (-1L, Typ_TmpGrp, 0);



  //----------------------------------------------------------------
  // Must display surface before boundary.
  if(!strchr(opts, 's') ) goto L_bvi;
    // display brkls->fia as shaded triangles
    GL_att_su (ATT_COL_GREEN);
    i2 = brkls->fiNr;
    ia = brkls->fia;
    for(i1=0; i1 < i2; ++i1) {
      // if(fa[ii].st < 0) continue;   // skip deleted faces
      UFA_Disp_fac2 (&MFAC__[ia[i1]], brkls->pta, zVal);
    }
    igl = 1;  // surfs



  //----------------------------------------------------------------
  L_bvi:
  if(!strchr(opts, 'b') ) goto L_glEnd;

    // boundary-right
    GL_att_cv (Typ_Att_hili1);   // INF_COL_CV
    igl = 0;  // pts, cvs, surfs

    for(i1=0; i1 < brkls->c1Nr; ++i1) {
      // get the points of edge cf1a[i1],ce1a[i1]
      UFA_2pt2_fac_esn (&pta[0], &pta[1],
                        brkls->ce1a[i1], &MFAC__[brkls->cf1a[i1]],
                        brkls->pta);
        // UT3D_stru_dump (Typ_PT2, &pta[0], " ps");
        // UT3D_stru_dump (Typ_PT2, &pta[1], " pe");
      GL_Disp_cv2z (2, pta, zVal);
    }


    // boundary-left
    GL_att_cv (Typ_Att_top1);

    for(i1 = brkls->c2Nr - 1; i1 >= 0; --i1) {
      // get the points of edge cf2a[i1],ce2a[i1]
      UFA_2pt2_fac_esn (&pta[0], &pta[1],
                        brkls->ce2a[i1], &MFAC__[brkls->cf2a[i1]],
                        brkls->pta);
      GL_Disp_cv2z (2, pta, zVal);
    }




  //----------------------------------------------------------------
  L_glEnd:
  GL_EndList1 (igl);     // glEndList


  //----------------------------------------------------------------
  L_end:


  return 0;

}


//==============================================================================
  int MSH2D_ctess__ (MemTab(Fac3) *fca, int *iba, int icNr, Point2 *p2a) {
//==============================================================================
// MSH2D_ctess__                  tesselate 2D-contour
// contour in is CCW, faces out are CCW
// Input:
//   p2a    2D-PunkteTabelle
//   iba    2D-Kontur (indices in die 2D-PunkteTabelle); not closed, CCW
//   icNr   nr of point-indices in iba
// Output:
//   fca    Faces (indices in die 2D-PunkteTabelle als *Fac3), CCW
//   retCod 0=OK;  -1=Err.
// 
// Fuer je 2 Kanten aufeinanderfolgende Kanten mit SpitzWinkel die Laenge
//   fuer Abkuerzung (dist p1-p3) merken.
//  - fuer gesamte Kontur.
//  Die beiden Kanten mit der kuerzesten Strecke suchen:
//    ein Face bilden,
//    die beiden Kanten durch eine Kante ersetzen (den Konturpunkt loeschen)
// 
// Testprog: MSH2D_test_ctess__ ();
// was MSH_itess_pt2n ../myAPPS/GIS__.c:4895


  int     i1, i2, ip0, ip1, ip2, ip3, ie1, ie2, ie3, iMin, *ica;
  double  d1, lMin;
  Point2  *pt1, *pt2;


  // printf("MSH2D_ctess__ %d\n",icNr);
  // for(i1=0;i1<icNr;++i1)UT3D_stru_dump (Typ_PT, &p2a[iba[i1]], "p-%d",i1);



  // remove last point if last-point == first-point
  pt1 = &p2a[iba[0]];
  pt2 = &p2a[iba[icNr-1]];
  if(UT2D_comp2pt(pt1,pt2,MSH2D_tol)) --icNr;
    // printf("  icNr=%d\n",icNr);


  //================================================================
  // OutsideContour 3 points only: single face. Done.
  if(icNr < 4) {
    if(icNr < 3) {TX_Print("MSH2D_ctess__ I001"); return -1;}
    // UFA_add_fac__ (MEMTAB_DAT(fca), ica[0], ica[1], ica[2]);
    UFA_add_fac__ (fca, iba[0], iba[1], iba[2]);
    goto L_exit;
  }


  //================================================================
  if(icNr >= 5) goto L_5_pts;
  L_4_pts:

  // 4 points - test convex
  i1 = MSH2D_ctess_4pt (iba, p2a);
  if(i1) {
    // -1 = faces 012 and 230
    UFA_add_fac__ (fca, iba[0], iba[1], iba[2]);
    UFA_add_fac__ (fca, iba[2], iba[3], iba[0]);
      // printf(" ctess__-fac: %d %d %d\n",iba[0], iba[1], iba[2]);
      // printf(" ctess__-fac: %d %d %d\n",iba[2], iba[3], iba[0]);
  } else {
    // 0 = faces 013 and 123
    UFA_add_fac__ (fca, iba[0], iba[1], iba[3]);
    UFA_add_fac__ (fca, iba[1], iba[2], iba[3]);
      // printf(" ctess__-fac: %d %d %d\n",iba[0], iba[1], iba[3]);
      // printf(" ctess__-fac: %d %d %d\n",iba[1], iba[2], iba[3]);
  }

    //TESTBLOCK
    // printf(" _itess_pt2n-4-points: i1=%d fNr=%d\n",i1,fca->rNr);
    // UFA_fac_dump1 (&fca->data[fca->rNr - 2], "%d",fca->rNr - 2);
    // UFA_fac_dump1 (&fca->data[fca->rNr - 1], "%d",fca->rNr - 1);
    // APT_disp_SymB (SYM_TRI_S, ATT_COL_BLUE, (Point2*)&p2a[ip2]);
    //TESTBLOCK

  goto L_exit;



  //================================================================
  L_5_pts:   // 5 or more points ..

  // copy iba -> ica
  i1 = icNr * sizeof(int);
  ica = (int*) MEM_alloc_tmp (i1);
  memcpy (ica, iba, i1);


  L_nxt_e5:
  ie1 = 0;
  lMin = UT_VAL_MAX;
  iMin = -1;
  iMin = 0;


  //----------------------------------------------------------------
  // loop tru boundary;
  L_nxt_e2:
  // ie1, ie2, ie3 are the next 3 points = 2 edges.
  ie2 = ie1 + 1;
  if(ie2 >= icNr) ie2 = 0;

  ie3 = ie2 + 1;
  if(ie3 >= icNr) ie3 = 0;

  ip1 = ica[ie1];
  ip2 = ica[ie2];
  ip3 = ica[ie3];

  // test dead end (can come up !)
  if(ie1 == ie3) {
    // printf("MSH2D_ctess__ I001 dead end\n");
    goto L_rem;
  }
    // printf(" ip1/2/3 = %d %d %d ie1/2/3 = %d %d %d icNr=%d\n",
           // ip1,ip2,ip3,ie1,ie2,ie3,icNr);


  // check if edges ie1-ie2-ie3 go left or right. contour is CCW.
  i1 = UT2D_sid_3pt ((Point2*)&p2a[ip1], (Point2*)&p2a[ip2], (Point2*)&p2a[ip3]);
    // printf(" ip1=%d ip2=%d ip3=%d side=%d\n",ip1,ip2,ip3,i1);
  // (i1 = 1): goes left; Ok to remove ip2
  // (i1 = -1)  goes right: continue (do not remove ip2)
  // (i1 = 0)   straight: continue (do not remove ip2)
  if(i1 <= 0) goto L_cont;  // L_nxt_e2;


  // test angle; select smaller angle 
  // d1 = UT2D_angr_3ptdr ((Point2*)&p2a[ip1],
                        // (Point2*)&p2a[ip3],
                        // (Point2*)&p2a[ip2], -1);
  // or before/behind
  // d1 = UT2D_sidPerp_3pt ((Point2*)&p2a[ip1],
                         // (Point2*)&p2a[ip2],
                         // (Point2*)&p2a[ip3]);
  // or parameter 
  UT2D_parpt_2pt (&d1, &p2a[ip1], &p2a[ip2], &p2a[ip3]);
    // printf(" d1=%f %d %d %d\n",d1,ip1,ip2,ip3);


  // select smaller value 
  if(d1 >= lMin) goto L_cont;
    // printf(" sel %d %d %d of %d\n",ip1,ip2,ip3,icNr);


  // test if other point falls into triangle ip1,ip2,ip3 
  // must be  left of seg2 AND seg3
  for(i1=0; i1<icNr; ++i1) {
    if((i1==ie1)||(i1==ie2)||(i1==ie3)) continue;
    i2 = UT2D_ck_pt_in_tria__ (&p2a[ip1], &p2a[ip2],
                               &p2a[ip3], &p2a[ica[i1]]);
      // printf(" p[%d] _in_ %d %d %d = %d\n",ica[i1],ip1,ip2,ip3,i2);
    if(i2 < 1) goto L_cont;
  }


  // select this point as to be removed if no better is found ..
  lMin = d1;
  iMin = ie1;


  L_cont:
  ++ie1;
  if(ie1 < icNr) goto L_nxt_e2;

  //----------------------------------------------------------------
  L_del_nxt:
  ie1 = iMin;

  ie2 = ie1 + 1;
  if(ie2 >= icNr) ie2 = 0;

  ie3 = ie2 + 1;
  if(ie3 >= icNr) ie3 = 0;
    // printf("+++ CreFac %d %d %d\n", ica[ie1], ica[ie2], ica[ie3]);
    // GR_Disp_ln2 (&p2a[ica[ie1]], &p2a[ica[ie3]], 8);

  // add face ie1-ie2-ie3.
    // GR_Disp_su3p2 (&p2a[ica[ie1]], &p2a[ica[ie2]], &p2a[ica[ie3]], 10);
  UFA_add_fac__ (fca, ica[ie1], ica[ie2], ica[ie3]);

    // TESTOUTPUT:
    // MSH_test_disp2_gc1 (fca->rNr - 1, fca->data, p2a);
    // printf("_add_ %d %d %d fNr=%d l=%lf\n",ica[ie1],ica[ie2],ica[ie3],
               // fca->rNr - 1,lMin);
    // printf(" del=%d l=%lf icNr=%d\n",ica[ie2],lMin,icNr);
    // if((fDbg)&&(icNr==23)) MSH__disp2_cv__ (ica, icNr, p2a, 8);


  // remove point ie2 from contour
  L_rem:
    // printf(" del pt %d\n",ie2);
  i1 = MEM_del_IndRec (&icNr, ica, ie2);
  if(i1 < 0) return -1;

    // TESTBLOCK
    // if(icNr < 9) goto L_exit; // TEST ONLY
    // END TESTBLOCK

  if(icNr > 4) goto L_nxt_e5;
  iba = ica;
  goto L_4_pts;


  //================================================================
  L_exit:
  // update recordNr of MSH2D_mtfnb
  // if(MSH2D_nfb_upd()) return -99;


    // TESTBLOCK
    // for(i1=0;i1<icNr;++i1)UT3D_stru_dump (Typ_PT, &p2a[ica[i1]], "p-%d",i1);
    // UFA_fac_dump_ (MEMTAB_DAT(fca), MEMTAB_IND(fca), "ex MSH2D_ctess__");
    // UFA_view__ (MEMTAB_DAT(fca), 0, MEMTAB_IND(fca), p2a, 0, "bf", -1L);
    // END TESTBLOCK

  return 0;

}


//================================================================
  int MSH2D_ctess_4pt (int *ia, Point2 *pa)  {
//================================================================
// check if indexed contour is convex | concave
// Output for pNr=4:
//   retCod:   0   faces 013 and 123
//            -1   faces 012 and 230
// see also UT2D_test_is_convex_inpt3

  int      i1, i2;



  // 4 points
  //          2                       1
  //         |  \                    |  \
  //         |    \                  |    \
  //         |      \                |      \
  //         |   0---1               |   3---0
  //         |  /                    |  /
  //         | /                     | /
  //          3                       2
  //       012 and 230            013 and 123

    // printf("MSH2D_ctess_4pt: %d %d %d %d\n",ia[0],ia[1],ia[2],ia[3]);

  // get i1 = if ia2 left or right of ia0-ia1
  i1 = UT2D_sid_3pt (&pa[ia[2]], &pa[ia[0]], (Point2*)&pa[ia[1]]);
  // get i2 = if ia0 left or right of ia2-ia3
  i2 = UT2D_sid_3pt (&pa[ia[0]], &pa[ia[2]], (Point2*)&pa[ia[3]]);
  // 0=parallel, 1=ip3_is_left_of_ip1-ip2, -1=right
    // printf(" _4_pts: i1=%d i2=%d\n",i1,i2);

  if((i1 > 0) && (i2 > 0)) return -1;  // faces 012 and 230
  else return 0;                       // faces 013 and 123

}

 
//================================================================
//================================================================
// Liste_TESTFUNKTIONEN:
//================================================================
//================================================================

#ifdef DEB

//================================================================
  int MSH2D_test_inp_CCW (TessStru *ts1) {
//================================================================
// test inputContours if CCW

  int      irc, i1, itNr, pNr, *ia, *ipa;
  double   d1;
  IndTab   *itTab, *itAct;
  Point2   *pa2;


  printf("MSH2D_test_inp_CCW \n");

  itNr  = MEMTAB_IND (&ts1->tab);  // nr of contours
  itTab = MEMTAB_DAT (&ts1->tab); // contours
  pa2 = MEMTAB_DAT (&ts1->pa2);
  ipa  = MEMTAB_DAT (&ts1->ipa);



  for(i1=0; i1<itNr; ++i1) {
    itAct = &itTab[i1];
    if((itAct->typi != MSH_EDGLN_OB) &&
       (itAct->typi != MSH_EDGLN_IB)) continue;
      // printf(" nxt-bnd i1=%d\n",i1);

    ia = &ipa[itAct->ibeg];       // first index
    pNr = itAct->iNr;             // nr of points
    // get sense-of-rotation; 1=CCW
    irc = UT2D_srar_inpt (&d1, pNr, ia, pa2);
      // printf(" srar_inpt %d %f\n",irc,d1);
    if(irc != 1) { TX_Error("Contour %d NOT CCW",i1); return -99; }
  }

  return 0;

}


//================================================================
  int MSH2D_test_ctess__ () {
//================================================================
// test tesselator for closed 2D-contours

#define NR_PTS 11

  int            *ica, i1;
  long           ld;
  MemTab(Point2) pa=_MEMTAB_NUL;
  MemTab(Fac3)   fa=_MEMTAB_NUL;
/*
  Point         OB1[] = {
  {0.514,    0.434,    0.},
  {0.474,    0.233,    0.},
  {0.589,    0.336,    0.},
  {0.773,    0.417,    0.},
  {0.66, 0.59, 0.},
  {0.60, 0.60, 0.},
  };
*/
  Point2        OB1[] = {
  { 0., 0.}, { 1., 0.}, { 2., 1.},
  { 2., 2.}, { 1., 3.}, { 0., 3.},
  {-1., 3.}, {-2., 2.}, {-2., 1.},
  {-1., 0.}, { 0., 0.}
  };


  // make fa = MemTab(Fac3)
  MemTab_ini (&fa, sizeof(Fac3), Typ_Fac3, NR_PTS * 4);

  // make p2a = table of points
  MemTab_ini (&pa, sizeof(Point), Typ_PT, NR_PTS);

  for(i1=0;i1<NR_PTS;++i1) GR_Disp_txi(&OB1[i1],i1,ATT_COL_RED);

  // add OB1 -> p2a
  MemTab_add (&pa, &ld, OB1, NR_PTS, 0);

  // make ica = table of indices of points
  ica = (int*) MEM_alloc_tmp (NR_PTS * (int)sizeof(int));
  for(i1=0;i1<NR_PTS;++i1) ica[i1] = i1;

  // tesselate ..
  MSH2D_tol = 0.005;
  MSH2D_ctess__ (&fa, ica, NR_PTS, OB1);

  // display
  UFA_view__ (MEMTAB_DAT(&fa),0,MEMTAB_IND(&fa),OB1,NR_PTS,"bf",-1L);

  return 0;
}

 
//================================================================
  int MSH2D_test__ () {
//================================================================

  int      irc, if1, if2, esn1, esn2;
  // Fac3     *fa;
  TessStru *ts1 = MSH2D_ts;

  // fa = MEMTAB_DAT (&ts1->fac);

  if1 = 0; if2 = 2;

/*
  // get common edges of 2 faces
  UFA_2esn_comm2fac (&esn1, &esn2, &fa[if1], &fa[if2]);
    printf(" if1=%d esn1=%d if2=%d esn2=%d\n",if1,esn1,if2,esn2);

  // check if esn1,esn2 is a common edge
  irc = UFA_2esn_comm_ck (esn1, esn2, &fa[if1], &fa[if2]);
  if(irc) return -1;
    printf(" irc=%d if1=%d esn1=%d if2=%d esn2=%d\n",irc,if1,esn1,if2,esn2);

  UFA_2fac_flip (if1, esn1, if2, esn2, fa, MSH2D_fnb);

  MSH2D_4fac_2fac_split (5, 1);
*/


  return 0;

}


//================================================================
  int MSH2D_test_qbx_disp (GridBox *qbx) {
//================================================================
// display gridlines

  int    i1, uNr, vNr;
  long   l1;
  double uMin, vMin, uMax, vMax, dx, dy, du, dv;
  Point2 p21, p22;


  uNr = qbx->ix;
  vNr = qbx->iy;

  du = qbx->dx;
  dv = qbx->dy;

  uMin = qbx->pMin.x;
  vMin = qbx->pMin.y;

  uMax = uMin + (uNr * du);
  vMax = vMin + (vNr * dv);

  printf("MSH2D_test_qbx_disp %f %f %f %f\n",uMin,vMin,uMax,vMax);
  printf(" uNr=%d vNr=%d\n",uNr,vNr);



  l1 = -1L;
  GL_DrawLn_Ini (&l1, 5);
  dx = uMin;
  dy = vMin;
  p21.y = vMin;
  p22.y = vMax;
  for(i1=0; i1<uNr; ++i1) {
    // next u = col = vert line
    p21.x = dx;
    p22.x = dx;
    GL_Disp_ln2 (&p21, &p22);
    dx += du;
  }
  p21.x = uMin;
  p22.x = uMax;
  for(i1=0; i1<vNr; ++i1) {
    // next row = hor line
    p21.y = dy;
    p22.y = dy;
    GL_Disp_ln2 (&p21, &p22);
    dy += dv;
  }
  GL_EndList ();

  return 0;

}


#endif


// EOF
