solmap_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.h"
#include "libmmg2d_private.h"
#include "mmg2dexterns_private.h"
#include "mmgexterns_private.h"

static inline
int MMG2D_solTruncatureForOptim(MMG5_pMesh mesh, MMG5_pSol met, int ani) {
  MMG5_pTria  ptt;
  MMG5_int    k;
  int         i,ier;
 
  assert ( mesh->info.optim || mesh->info.hsiz > 0. );
 
  /* Detect the points not used by triangles */
  ++mesh->base;
#ifndef NDEBUG
  for (k=1; k<=mesh->np; k++) {
    assert ( mesh->point[k].flag < mesh->base );
  }
#endif
 
  for (k=1; k<=mesh->nt; k++) {
    ptt = &mesh->tria[k];
    if ( !MG_EOK(ptt) ) continue;
 
    for (i=0; i<3; i++) {
      mesh->point[ptt->v[i]].flag = mesh->base;
    }
  }
 
  /* Compute hmin/hmax on unflagged points and truncate the metric */
  if ( !ani ) {
    ier = MMG5_solTruncature_iso(mesh,met);
  }
  else {
    MMG5_solTruncature_ani = MMG5_2dSolTruncature_ani;
    ier = MMG5_solTruncature_ani(mesh,met);
  }
 
  return ier;
}

int MMG2D_doSol_iso(MMG5_pMesh mesh,MMG5_pSol sol) {
  MMG5_pTria      ptt,pt;
  MMG5_pPoint     p1,p2;
  double          ux,uy,dd;
  MMG5_int        k,ipa,ipb;
  int             i,ib;
  int             MMG_inxtt[5] = {0,1,2,0,1};
 
  // here we guess that we have less than int32 edges passing through a point
  int             *mark;
  MMG5_SAFE_CALLOC(mark,mesh->np+1,int,return 0);
 
  /* Memory alloc */
  if ( sol->size!=1 ) {
    fprintf(stderr,"\n  ## Error: %s: unexpected size of metric: %d.\n",
            __func__,sol->size);
    return 0;
  }
 
  if ( !MMG2D_Set_solSize(mesh,sol,MMG5_Vertex,mesh->np,sol->size) )
    return 0;
 
  /* Travel the triangles edges and add the edge contribution to edges
   * extermities */
  for (k=1; k<=mesh->nt; k++) {
    ptt = &mesh->tria[k];
    if ( !ptt->v[0] )  continue;
 
    for (i=0; i<3; i++) {
      ib  = MMG_inxtt[i+1];
      ipa = ptt->v[i];
      ipb = ptt->v[ib];
      p1  = &mesh->point[ipa];
      p2  = &mesh->point[ipb];
 
      ux  = p1->c[0] - p2->c[0];
      uy  = p1->c[1] - p2->c[1];
      dd  = sqrt(ux*ux + uy*uy);
 
      sol->m[ipa] += dd;
      mark[ipa]++;
      sol->m[ipb] += dd;
      mark[ipb]++;
    }
  }
 
  /* vertex size */
  for (k=1; k<=mesh->np; k++) {
    if ( !mark[k] )  {
      continue;
    }
    sol->m[k] = sol->m[k] / (double)mark[k];
  }
  MMG5_SAFE_FREE(mark);
 
  /* Computation of hmin/hmax if not provided + size truncature */
  MMG2D_solTruncatureForOptim(mesh,sol,0);
 
  /* compute quality */
  if ( MMG2D_caltri ) {
    for (k=1; k<=mesh->nt; k++) {
      pt = &mesh->tria[k];
      pt->qual = MMG2D_caltri_iso(mesh,sol,pt);
    }
  }
 
  return 1;
}

int MMG2D_doSol_ani(MMG5_pMesh mesh,MMG5_pSol sol) {
  MMG5_pTria      ptt,pt;
  MMG5_pPoint     p1,p2;
  double          ux,uy,dd,tensordot[3];
  MMG5_int        k,ipa,ipb,iadr;
  int             i,ib;
  int             MMG_inxtt[5] = {0,1,2,0,1};
 
  // here we guess that we have less than int32 edges passing through a point
  int             *mark;
  MMG5_SAFE_CALLOC(mark,mesh->np+1,int,return 0);
 
  /* Memory alloc */
  if ( sol->size!=3 ) {
    fprintf(stderr,"\n  ## Error: %s: unexpected size of metric: %d.\n",
              __func__,sol->size);
      return 0;
    }
 
  if ( !MMG2D_Set_solSize(mesh,sol,MMG5_Vertex,mesh->np,sol->size) )
    return 0;
 
  /* Travel the triangles edges and add the edge contribution to edges
   * extermities */
  for (k=1; k<=mesh->nt; k++) {
    ptt = &mesh->tria[k];
    if ( !ptt->v[0] )  continue;
 
    for (i=0; i<3; i++) {
      ib  = MMG_inxtt[i+1];
      ipa = ptt->v[i];
      ipb = ptt->v[ib];
      p1  = &mesh->point[ipa];
      p2  = &mesh->point[ipb];
 
      ux  = p1->c[0] - p2->c[0];
      uy  = p1->c[1] - p2->c[1];
 
      tensordot[0] = ux*ux;
      tensordot[1] = ux*uy;
      tensordot[2] = uy*uy;
 
      iadr = 3*ipa;
      sol->m[iadr]   += tensordot[0];
      sol->m[iadr+1] += tensordot[1];
      sol->m[iadr+2] += tensordot[2];
      mark[ipa]++;
 
      iadr = 3*ipb;
      sol->m[iadr]   += tensordot[0];
      sol->m[iadr+1] += tensordot[1];
      sol->m[iadr+2] += tensordot[2];
      mark[ipb]++;
    }
  }
 
  /* Compute metric tensor and hmax if not specified */
  for (k=1; k<=mesh->np; k++) {
    if ( !mark[k] ) {
      continue;
    }
 
    iadr = 3*k;
    /* Metric = nedges/dim * inv (sum(tensor_dot(edges,edges))).
     * sum(tensor_dot) is stored in sol->m so reuse tensordot to
     * compute M.  */
    dd = 1./(sol->m[iadr]*sol->m[iadr+2] - sol->m[iadr+1]*sol->m[iadr+1]);
    dd *= (double)mark[k]*0.5;
 
    tensordot[0] = sol->m[iadr+2];
    tensordot[1] = -sol->m[iadr+1];
    tensordot[2] = sol->m[iadr];
 
    sol->m[iadr]   = dd*tensordot[0];
    sol->m[iadr+1] = dd*tensordot[1];
    sol->m[iadr+2] = dd*tensordot[2];
 
#ifndef NDEBUG
    /* Check metric */
    double lambda[2],vp[2][2];
    MMG5_eigensym(sol->m+iadr,lambda,vp);
 
    assert (lambda[0] > 0. && lambda[1] > 0. && "Negative eigenvalue");
 
    /* Normally the case where the point belongs to only 2 colinear points is
    impossible */
    assert (isfinite(lambda[0]) && isfinite(lambda[1]) && "wrong eigenvalue");
#endif
  }
  MMG5_SAFE_FREE(mark);
 
  /* Size truncature */
  MMG2D_solTruncatureForOptim(mesh,sol,1);
 
  /* compute quality */
  if ( MMG2D_caltri ) {
    for (k=1; k<=mesh->nt; k++) {
      pt = &mesh->tria[k];
      pt->qual = MMG2D_caltri_ani(mesh,sol,pt);
    }
  }
 
  return 1;
}