isosiz_2d.c

Go to the documentation of this file.

/* =============================================================================
**  This file is part of the mmg software package for the tetrahedral
**  mesh modification.
**  Copyright (c) Bx INP/CNRS/Inria/UBordeaux/UPMC, 2004-
**
**  mmg is free software: you can redistribute it and/or modify it
**  under the terms of the GNU Lesser General Public License as published
**  by the Free Software Foundation, either version 3 of the License, or
**  (at your option) any later version.
**
**  mmg is distributed in the hope that it will be useful, but WITHOUT
**  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
**  FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
**  License for more details.
**
**  You should have received a copy of the GNU Lesser General Public
**  License and of the GNU General Public License along with mmg (in
**  files COPYING.LESSER and COPYING). If not, see
**  <http://www.gnu.org/licenses/>. Please read their terms carefully and
**  use this copy of the mmg distribution only if you accept them.
** =============================================================================
*/
#include "libmmg2d_private.h"
#include "libmmg2d.h"
#include "mmg2dexterns_private.h"
#include "mmgexterns_private.h"

int MMG2D_sum_reqEdgeLengthsAtPoint(MMG5_pMesh mesh,MMG5_pSol met,MMG5_pTria pt,int8_t i) {
  MMG5_int         ip0,ip1;
 
  ip0 = pt->v[MMG5_iprv2[i]];
  ip1 = pt->v[MMG5_inxt2[i]];
 
  if ( !MMG5_sum_reqEdgeLengthsAtPoint(mesh,met,ip0,ip1) )
    return 0;
 
  return 1;
}
 

int MMG2D_set_metricAtPointsOnReqEdges ( MMG5_pMesh mesh,MMG5_pSol met, int8_t ismet ) {
  MMG5_pTria pt;
  MMG5_int   k,iadj;
  int        i;
 
  /* Reset the tria flag */
  for ( k=1; k<=mesh->nt; k++ ) {
    mesh->tria[k].flag = 0;
  }
 
  /* Reset the input metric at required edges extremities */
  if ( !MMG5_reset_metricAtReqEdges_surf (mesh, met, ismet ) ) {
    return 0;
  }
 
  /* Process the required edges and add the edge length to the metric of the
   * edge extremities */
  for ( k=1; k<=mesh->nt; k++ ) {
    pt = &mesh->tria[k];
    if ( !MG_EOK(pt) )  continue;
 
    /* Mark the tria as proceeded */
    pt->flag = 1;
    for ( i=0; i<3; i++ ) {
      if ( (pt->tag[i] & MG_REQ) || (pt->tag[i] & MG_NOSURF) ||
           (pt->tag[i] & MG_PARBDY) ) {
        /* Check if the edge has been proceeded by the neighbour triangle */
        iadj = mesh->adja[3*(k-1)+i+1];
        if ( iadj && mesh->tria[iadj/3].flag ) continue;
 
        if ( !MMG2D_sum_reqEdgeLengthsAtPoint(mesh,met,pt,i) ) {
          return 0;
        }
      }
    }
  }
 
  /* Travel the points and compute the metric of the points belonging to
   * required edges as the mean of the required edges length */
  if ( !MMG5_compute_meanMetricAtMarkedPoints ( mesh,met ) ) {
    return 0;
  }
 
return 1;
}

int MMG2D_defsiz_iso(MMG5_pMesh mesh,MMG5_pSol met) {
  MMG5_pTria  pt;
  MMG5_pPoint p0,p1,p2;
  MMG5_pPar   ppa;
  double      t1[2],t2[2],b1[2],b2[2],gpp1[2],gpp2[2],pv,M1,M2;
  double      ps1,ps2,ux,uy,ll,li,lm,hmax,hausd,hmin,lhmax,lhausd;
  MMG5_int    k,ip,ip1,ip2;
  int         l;
  int8_t      ismet;
  uint8_t     i,i1,i2;
 
  if ( !MMG5_defsiz_startingMessage (mesh,met,__func__) ) {
    return 0;
  }
 
  /* Reset the s and flag field : s is used to count the number of req edges to
   * which a req point belongs and flag is used to differenciate the point not
   * yet treated, the points with a required metric and the classical points */
  for (k=1; k<=mesh->np; k++) {
    p0 = &mesh->point[k];
    p0->flag = 0;
    p0->s    = 0;
  }

  hmax = mesh->info.hmax;
  hausd = mesh->info.hausd;
  hmin = mesh->info.hmin;
 
  /* Allocate the structure */
  if ( !met->np ) {
    ismet = 0;
 
    /* Allocate and store the header informations for each solution */
    if ( !MMG2D_Set_solSize(mesh,met,MMG5_Vertex,mesh->np,MMG5_Scalar) ) {
      return 0;
    }
  }
  else {
    ismet = 1;
    assert ( met->m );
  }

  if ( !mesh->info.nosizreq ) {
    if ( !MMG2D_set_metricAtPointsOnReqEdges ( mesh,met,ismet ) ) {
      return 0;
    }
  }

  if ( !ismet ) {
    /* Initialize metric with a constant size */
    for ( k=1; k<=mesh->np; k++ ) {
      if ( mesh->point[k].flag ) {
        continue;
      }
      met->m[k] = hmax;
      mesh->point[k].flag = 1;
    }
  }

  for ( k=1; k<=mesh->nt; k++ ) {
    pt = &mesh->tria[k];
    if ( !MG_EOK(pt) ) continue;
 
    for ( i=0; i<3; i++ ) {
 
      if ( !MG_EDG(pt->tag[i]) ) continue;
 
      i1 = MMG5_inxt2[i];
      i2 = MMG5_iprv2[i];
      ip1 = pt->v[i1];
      ip2 = pt->v[i2];
 
      p1 = &mesh->point[ip1];
      p2 = &mesh->point[ip2];
 
      if ( p1->flag>1 && p2->flag>1 ) continue;
 
      lhmax = hmax;
      lhausd = hausd;
 
      /* Retrieve local parameters associated to edge i */
      if ( mesh->info.npar ) {
        for (l=0; l<mesh->info.npar; l++) {
          ppa = &mesh->info.par[l];
          if ( ppa->elt == MMG5_Edg && ppa->ref == pt->edg[i] ) {
            lhmax = ppa->hmax;
            lhausd = ppa->hausd;
            break;
          }
        }
      }
 
 
      ux = p2->c[0] - p1->c[0];
      uy = p2->c[1] - p1->c[1];
      ll = ux*ux + uy*uy;
      if ( ll < MMG5_EPSD ) continue;
      li = 1.0 / sqrt(ll);
 
      /* Recovery of the two tangent vectors associated to points p1,p2; they
       * need not be oriented in the same fashion */
      if ( MG_CRN & p1->tag || (p1->tag & MG_NOM) ) {
        t1[0] = li*ux;
        t1[1] = li*uy;
      }
      else {
        t1[0] = -p1->n[1];
        t1[1] = p1->n[0];
      }
 
      if ( MG_CRN & p2->tag || (p2->tag & MG_NOM) ) {
        li = 1.0 / sqrt(ll);
        t2[0] = li*ux;
        t2[1] = li*uy;
      }
      else {
        t2[0] = -p2->n[1];
        t2[1] = p2->n[0];
      }
 
      /* Calculation of the two Bezier coefficients along the curve */
      ps1   = ux*t1[0] + uy*t1[1];
      b1[0] = p1->c[0] + MMG5_ATHIRD*ps1*t1[0];
      b1[1] = p1->c[1] + MMG5_ATHIRD*ps1*t1[1];
 
      ps2   = ux*t2[0]+uy*t2[1];
      b2[0] = p2->c[0] - MMG5_ATHIRD*ps2*t2[0];
      b2[1] = p2->c[1] - MMG5_ATHIRD*ps2*t2[1];
 
      ps1 *= ps1;
      ps2 *= ps2;
 
      if ( ps1 < MMG5_EPSD || ps2 < MMG5_EPSD ) continue;
 
      /* \gamma^{\prime\prime}(0); \gamma^\prime(0) = ps*t1 by construction */
      gpp1[0] = 6.0*(p1->c[0] - 2.0*b1[0] + b2[0]);
      gpp1[1] = 6.0*(p1->c[1] - 2.0*b1[1] + b2[1]);
 
      /* Vector product gpp1 ^ t1 */
      pv = gpp1[0]*t1[1] - gpp1[1]*t1[0];
      M1 = fabs(pv)/ps1;
 
      /* \gamma^{\prime\prime}(1); \gamma^\prime(1) = -ps*t2 by construction */
      gpp2[0] = 6.0*(p2->c[0] - 2.0*b2[0] + b1[0]);
      gpp2[1] = 6.0*(p2->c[1] - 2.0*b2[1] + b1[1]);
 
      /* Vector product gpp2 ^ t2 */
      pv = gpp2[0]*t2[1] - gpp2[1]*t2[0];
      M2 = fabs(pv)/ps2;
 
      M1 = MG_MAX(M1,M2);
      if ( M1 < MMG5_EPSD )
        lm = lhmax;
      else {
        lm = 8.0*lhausd / M1;
        lm = MG_MIN(lhmax,sqrt(lm));
      }
 
      if ( p1->flag < 3 ) {
        met->m[ip1] = MG_MAX(hmin,MG_MIN(met->m[ip1],lm));
      }
      if ( p2->flag < 3 ) {
        met->m[ip2] = MG_MAX(hmin,MG_MIN(met->m[ip2],lm));
      }
    }
  }

  /* Without local parameters information, only the boundary edges impose a
   * minimum size feature */
  if ( mesh->info.npar ) {
    /* Minimum size feature imposed by triangles */
    for (k=1; k<=mesh->nt; k++) {
      pt = &mesh->tria[k];
      if ( !MG_EOK(pt) ) continue;
 
      /* Retrieve local parameters associated to triangle k */
      for (l=0; l<mesh->info.npar; l++) {
        ppa = &mesh->info.par[l];
        if ( ppa->elt == MMG5_Triangle && ppa->ref == pt->ref ) {
          for (i=0; i<3; i++) {
            ip = pt->v[i];
            if ( mesh->point[ip].flag < 3 ) {
              met->m[ip] = MG_MAX(hmin,MG_MIN(met->m[ip],ppa->hmax));
            }
          }
          break;
        }
      }
    }
    /* Minimum size feature imposed by vertices */
    for (k=1; k<=mesh->np; k++) {
      p0 = &mesh->point[k];
      if ( (!MG_VOK(p0)) || p0->flag == 3 ) continue;
 
      /* Retrieve local parameters associated to vertex k */
      for (l=0; l<mesh->info.npar; l++) {
        ppa = &mesh->info.par[l];
        if ( ppa->elt == MMG5_Vertex && ppa->ref == p0->ref ) {
          met->m[k] = MG_MAX(hmin,MG_MIN(met->m[k],ppa->hmax));
          break;
        }
      }
    }
  }
 
  return 1;
}