Bamgx.cpp

Go to the documentation of this file.

 
#include "./Bamgx.h"
#include "../../bamg/bamgobjects.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
 
using namespace bamg;
using namespace std;
 
int Bamgx(BamgMesh* bamgmesh_out,BamgGeom* bamggeom_out,BamgMesh* bamgmesh_in,BamgGeom* bamggeom_in,BamgOpts* bamgopts){
 
   /*Bamg options*/
   int    maxnbv;
   double coef;
   int    verbosity;
   int    nbsmooth;
 
   /*intermediary*/
   int i;
   int noerr=1;
   double hminaniso=1e-100;
   Mesh* Thr=NULL;
   Mesh* Thb=NULL;
 
   /*Bamg options*/
   nbsmooth =bamgopts->nbsmooth;
   coef     =bamgopts->coeff;
   maxnbv   =bamgopts->maxnbv;
   verbosity=bamgopts->verbose;
 
   // no metric -> no smoothing
   if (bamgopts->metric==NULL) nbsmooth=0;
 
   /*If no mesh in input, generate one*/
   if(bamgmesh_in->TrianglesSize[0]==0){
      /*Mesh generation {{{*/
 
      //Step1: generate geometry Gh
      if (verbosity>0) _printf_("Construction of a mesh from a given geometry\n");
      if (verbosity>1) _printf_("   Processing geometry...\n");
      Geometry Gh(bamggeom_in,bamgopts);
 
      //get hmin and hmax from geometry to generate the metric
      bamgopts->hmin = Max(bamgopts->hmin,Gh.MinimalHmin());
      bamgopts->hmax = Min(bamgopts->hmax,Gh.MaximalHmax());
 
      //build metric using geometry
      if (verbosity>1) _printf_("   Generating Metric...\n");
      for(i=0;i<Gh.nbv;i++){
         Metric M=Gh[i];
         EigenMetric Vp(M/coef);
         Vp.Maxh(bamgopts->hmax);
         Vp.Minh(bamgopts->hmin);
         Gh.vertices[i].m = Vp;
      }
 
      //generate mesh
      if (verbosity>1) _printf_("   Generating Mesh...\n");
      Mesh Th(maxnbv,Gh,bamgopts);
 
      //Split corners if requested
      if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();
 
      //Renumbering
      Th.TrianglesRenumberBySubDomain();
 
      //Crack mesh if requested
      if(bamgopts->Crack) Th.CrackMesh(bamgopts);
 
      //Build output
      if (verbosity>1) _printf_("   Write Mesh...\n");
      Th.WriteMesh(bamgmesh_out,bamgopts);
      if (verbosity>1) _printf_("   Write Geometry...\n");
      Gh.WriteGeometry(bamggeom_out,bamgopts);
 
      //clean up
   // delete &Th;
   // delete &Gh;
      /*}}}*/
   }
   else{
      /*Anisotropic mesh adaptation {{{*/
 
      // read background mesh
      if (verbosity>0) _printf_("Anisotropic mesh adaptation\n");
      if (verbosity>1) _printf_("   Processing initial mesh and geometry...\n");
      Mesh BTh(bamggeom_in,bamgmesh_in,bamgopts);
 
      //Make Quadtree from background mesh
      BTh.MakeBamgQuadtree();
 
      //Bound hmin and hmax
      bamgopts->hmin=Max(bamgopts->hmin,BTh.MinimalHmin());
      bamgopts->hmax=Min(bamgopts->hmax,BTh.MaximalHmax());
 
      //Generate initial metric
      if (bamgopts->metric){
         if (verbosity>1) _printf_("   Processing Metric...\n");
         BTh.ReadMetric(bamgopts);
      }
      else {
         if (verbosity>1) _printf_("   Generating initial Metric...\n");
         Metric Mhmax(bamgopts->hmax);
         for (int iv=0;iv<BTh.nbv;iv++) BTh[iv].m = Mhmax;
      }
 
      //use present fields to generate metric if present
      if (bamgopts->field){
         if (verbosity>1) _printf_("   Merge metric with field provided...\n");
         BTh.AddMetric(bamgopts);
      }
 
      // change using hVertices if provided
      if(bamgopts->hVertices && bamgopts->hVerticesLength==BTh.nbv){
         if (verbosity>1) _printf_("   Merging Metric with hVertices...\n");
         for (i=0;i<BTh.nbv;i++){
            if (!xIsNan<IssmPDouble>(bamgopts->hVertices[i])){
               BTh[i].m=Metric((float)bamgopts->hVertices[i]);
            }
         }
      }
 
      // change using hminVertices if provided
      if (bamgopts->hminVertices){
         if (verbosity>1) _printf_("   Merging Metric with hminVertices...\n");
         for (i=0;i<BTh.nbv;i++){
            if (!xIsNan<IssmPDouble>(bamgopts->hminVertices[i])){
               Metric M=BTh.vertices[i].m;
               EigenMetric Vp(M/coef);
               Vp.Minh(bamgopts->hminVertices[i]);
               BTh.vertices[i].m=Vp;
            }
         }
      }
 
      // change using hmaxVertices if provided
      if (bamgopts->hmaxVertices){
         if (verbosity>1) _printf_("   Merging Metric with hmaxVertices...\n");
         for (i=0;i<BTh.nbv;i++){
            if (!xIsNan<IssmPDouble>(bamgopts->hmaxVertices[i])){
               Metric M=BTh.vertices[i].m;
               EigenMetric Vp(M/coef);
               Vp.Maxh(bamgopts->hmaxVertices[i]);
               BTh.vertices[i].m=Vp;
            }
         }
      }
 
      //Smoothe metric
      BTh.SmoothMetric(bamgopts,bamgopts->gradation);
 
      //Control element subdivision
      BTh.MaxSubDivision(bamgopts,bamgopts->maxsubdiv);
 
      //Bound anisotropy
      BTh.BoundAnisotropy(bamgopts,bamgopts->anisomax,hminaniso);
 
      //Build new mesh
      if (verbosity>1) _printf_("   Generating Mesh...\n");
      Thr=&BTh,Thb=0;
      Mesh & Th( *(0 ?  new Mesh(*Thr,&Thr->Gh,Thb,maxnbv) :  new Mesh(maxnbv,BTh,bamgopts,bamgopts->KeepVertices)));
      if (Thr != &BTh) delete Thr;
 
      //Split corners if requested
      if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();
 
      //Renumber by subdomain
      Th.TrianglesRenumberBySubDomain();
 
      //Smooth vertices
      if(nbsmooth>0) Th.SmoothingVertex(bamgopts,nbsmooth,bamgopts->omega);
 
      //display info
      if(verbosity>0) {
         if (Th.nbt-Th.nbtout){
            _printf_("   new number of triangles = " << (Th.nbt-Th.nbtout) << "\n");
         }
      }
 
      //Build output
      if (verbosity>1) _printf_("   Write Mesh...\n");
      Th.WriteMesh(bamgmesh_out,bamgopts);
      if (verbosity>1) _printf_("   Write Geometry...\n");
      Th.Gh.WriteGeometry(bamggeom_out,bamgopts);
      if (verbosity>1) _printf_("   Write Metric...\n");
      BTh.WriteMetric(bamgopts);
 
      /*clean up*/
      delete &Th;
      //delete &BTh;
      /*}}}*/
   }
 
   /*No error return*/
   if (verbosity>1) _printf_("   Exiting Bamg.\n");
   return noerr;
 
}