bamg.cpp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

// ********** DO NOT REMOVE THIS BANNER **********
//
// SUMMARY:  Bamg: Bidimensional Anisotrope Mesh Generator
// RELEASE: 0 
//
// AUTHOR:   F. Hecht,    
// ORG    :  UMPC
// E-MAIL :   [e-mail]   
//
// ORIG-DATE:     Dec 97
// Modif          March 98
/*
 
 This file is part of Bamg or FreeFEm++
 
 Freefem++ is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.
 
 Freefem++  is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */
// E-MAIL :   [e-mail]   
//
// ORIG-DATE:     Dec 97
#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <setjmp.h>
#include <new>
#include <cassert>
#include "Meshio.h"
#include <iomanip>
#include "Mesh2.h"
#include "QuadTree.h"
using namespace bamg;
using namespace std;
#include <fstream>
#ifdef __MWERKS__
#define   NBVMAX 10000
#else
#define   NBVMAX 50000   
#endif 
//                       0         1         2         3
//                       0123456789012345678901234567890
const char whatbamg []= "@(#)LJLL, le 20 fevr 2006,  Bamg v1.02  ";
const char *bamgversion = whatbamg + 28;
int initgraph=0;
//long verbosity=2;

#ifdef __MWERKS__
#define R_OK 0
#define F_OK 0
#define W_OK 0
#include <unix.h>
inline int access( char *fileName, int notUsed )
{
	struct stat		statRec;
	return  stat( fileName, &statRec );
}
#endif
#ifdef _MSC_VER
#define R_OK 04
#define F_OK 0
#define W_OK 02
#include <io.h>
inline int access( char *fileName, int mode)
{
  return  _access( fileName, mode );
}
#endif

void NewHandler();
void NewHandler()
{ cerr << " Not enought memory " << endl;
  exit(1);
}
void forDebug();


void  MeshErrorIO(ios& )
{
  MeshError(999);
  exit(1);
}

#ifdef DRAWING 
bool withrgraphique= initgraphique; 
#else
bool withrgraphique=false; 
#endif

void forDebug()
{ 
#ifdef DRAWING 
if(initgraph) {
  if (CurrentTh)     CurrentTh->inquire();
  closegraphique();initgraph=0; }
#endif
}

int main(int argc, char **argv)
{  
#ifdef __MWERKS__
  cout << " version " << bamgversion << " for initialisation of the toolbox " << endl;
#endif 
  // JUST for HP CC bug 
#ifdef LONGLONG
  long long bidon=2;
  if( (double) bidon != 2.0) {
    cerr << "Compiler bug  Pb of cast of long long to a double => remove the LONGLONG define" << endl;
    exit(10); }
#endif
  double time0, time1,time2,time3,time4;
  //   2 way for uses ---
  //   1 first mesh 
  //   or adaptation 
  MeshIstreamErrorHandler = MeshErrorIO;
  //  double *sol=0;
  //  int adaptation = 0 ;
  Int4 i;
  hinterpole=1; // by def interpolation a h 
  //  int Err=0;
  int fileout=0;
  int nbvx = NBVMAX;
  int iso =0,AbsError=0,nbjacoby=1,allquad=0;
  int NoMeshReconstruction=0;
  int Rescaling = 1;
  double costheta=2;
  Real8 cutoffradian=-1;
  double anisomax = 1e6;
  double err=0.01,errg=0.1,coef=1,hmin=1.e-100,hmax=1.e17,raison=0,CutOff=1e-5;
  int KeepBackVertices=1;
 double hminaniso=1e-100; 
  const double  boundmaxsubdiv = 10;
  double  maxsubdiv=boundmaxsubdiv;
  double omega=1.8;
  int NbSmooth=3;
  double *solMbb =0,*solMBB=0;
  int  *typesolsBB =0;
  Int4 nbsolbb=0,lsolbb=0;
  Int4 nbsolBB=0,lsolBB=0;
  int SplitEdgeWith2Boundary=0;
  int rbbeqMbb=0,rBBeqMBB=0;
  int ChoiseHessien = 0;
  double power=1;
  Triangles *Thr=0,*Thb=0;

  char * fgeom=0,*fmeshback=0,*fmeshout=0,*fmeshr=0,*fmetrix=0,*famfmt=0,*fmsh=0,
    *fnopo=0,*fftq=0,*fam=0,*famdba=0,*rbb=0,*rBB=0,*wbb=0,*wBB=0,
    *fMbb=0,*foM=0,*fMBB=0;
  verbosity =2;
  char *datargv[128] ;
  int datargc=1;
  datargv[0]= datargv[1]=0;// for create a error if no parameter 
  const char * datafile = argc ==2 ? argv[1] : "DATA_bamg";
  
  atexit( forDebug);
  set_new_handler( &NewHandler);

  if (argc<3)  
    { // for test on the mac
      if (!access(datafile,R_OK))
	{
	  MeshIstream data(datafile);
	  while(data.cm() && datargc <128)
	    datargv[datargc++]=data.ReadStr();
	  
      if (argc==1)
	datargv[0]=argv[0]; // copy of the program name 
      argc=datargc;
      argv=datargv;
	}
    };

  for (i=1;i<argc;i++)
    if (!strcmp(argv[i],"-v")  && ++i<argc)
      verbosity = atoi(argv[i]);
    else if (!strcmp(argv[i],"-NoRescaling") ) 
       Rescaling = 0;
    else if (!strcmp(argv[i],"-Rescaling") ) 
       Rescaling = 1;
    else if (!strcmp(argv[i],"-H")  && ++i<argc) 
       ChoiseHessien = atoi(argv[i]);
    else if (!strcmp(argv[i],"-power")  && ++i<argc) 
       power = atof(argv[i]);
    else if (!strcmp(argv[i],"-nbv")   && ++i<argc)
      nbvx = atoi(argv[i]);
    else if (!strcmp(argv[i],"-nbs")   && ++i<argc)
       nbvx = atoi(argv[i]);
    else if (!strcmp(argv[i],"-g")   && ++i<argc)
      fgeom = argv[i];
    else if ((!strcmp(argv[i],"-thetaquad") || (!strcmp(argv[i],"-ThetaQuad"))) && ++i<argc)
      costheta  = 1-Abs(cos(Pi*atof(argv[i])/180.0));
    else if (!strcmp(argv[i],"-b")  && ++i<argc)
      fmeshback = argv[i];
    else if (!strcmp(argv[i],"-r")  && ++i<argc)
      fmeshr = argv[i];
    else if (!strcmp(argv[i],"-M")  && ++i<argc) 
      fmetrix = argv[i];
    else if (!strcmp(argv[i],"-rbb") && ++i<argc) 
      rbb = argv[i];
    else if (!strcmp(argv[i],"-rBB") && ++i<argc) 
      rBB = argv[i];
    else if (!strcmp(argv[i],"-Mbb") && ++i<argc) 
      fMbb = argv[i];
    else if (!strcmp(argv[i],"-MBB") && ++i<argc) 
      fMBB = argv[i];
    else if (!strcmp(argv[i],"-wBB") && ++i<argc) 
      wBB = argv[i];
    else if (!strcmp(argv[i],"-wbb") && ++i<argc) 
      wbb = argv[i];
     else if (!strcmp(argv[i],"-o")  && ++i<argc) 
       fmeshout= argv[i];
    else if (!strcmp(argv[i],"-intm")  && i<argc)
      hinterpole=0; 
    else if (!strcmp(argv[i],"-oam_fmt") && ++i<argc)
      famfmt = argv[i];
    else if (!strcmp(argv[i],"-oam") && ++i<argc)
      fam= argv[i];
    else if (!strcmp(argv[i],"-omsh") && ++i<argc)
      fmsh = argv[i];
    else if (!strcmp(argv[i],"-oftq") && ++i<argc)
      fftq = argv[i];
    else if (!strcmp(argv[i],"-onopo") && ++i<argc)
      fnopo = argv[i];
    else if (!strcmp(argv[i],"-oamdba") && ++i<argc)
      famdba= argv[i];
    else if (!strcmp(argv[i],"-oM") && ++i<argc)
      foM = argv[i];
    else if (!strcmp(argv[i],"-coef")  && ++i<argc) 
      coef = atof(argv[i]);
    else if (!strcmp(argv[i],"-err")  && ++i<argc) 
      err = atof(argv[i]);
    else if (!strcmp(argv[i],"-errg")  && ++i<argc) 
      errg = atof(argv[i]);
    else if (!strcmp(argv[i],"-raison")  && ++i<argc) 
      raison = atof(argv[i]);
    else if (!strcmp(argv[i],"-ratio")  && ++i<argc) 
      raison = atof(argv[i]);
    else if (!strcmp(argv[i],"-hmin")  && ++i<argc) 
      hmin= atof(argv[i]);
     else if (!strcmp(argv[i],"-hminaniso")  && ++i<argc) 
      hminaniso= atof(argv[i]);
   else if (!strcmp(argv[i],"-hmax")  && ++i<argc) 
       hmax = atof(argv[i]);
    else if (!strcmp(argv[i],"-thetamax")  && ++i<argc) 
      cutoffradian = atof(argv[i])*Pi/180.0;
    else if (!strcmp(argv[i],"-anisomax")  && ++i<argc) 
       anisomax =  atof(argv[i]);
    else if (!strcmp(argv[i],"-maxsubdiv")  && ++i<argc) 
       maxsubdiv = atof(argv[i]);
    else if (!strcmp(argv[i],"-iso"))
      iso = 1;
    else if (!strcmp(argv[i],"-KeepBackVertices"))
      KeepBackVertices = 1;
    else if (!strcmp(argv[i],"-noKeepBackVertices"))
      KeepBackVertices = 0;
    else if (!strcmp(argv[i],"-KeepBackVertices"))
	KeepBackVertices = 1;
    else if (!strcmp(argv[i],"-noKeepBackVertices"))
	KeepBackVertices = 0;
    else if (!strcmp(argv[i],"-2q") || !strcmp(argv[i],"-2Q"))
      allquad=1;
    else if (!strcmp(argv[i],"-2") )
      allquad=2;
    else if (!strcmp(argv[i],"-aniso"))
      iso = 0;
    else if (!strcmp(argv[i],"-RelError"))
      AbsError = 0;
    else if (!strcmp(argv[i],"-AbsError"))
      AbsError = 1;
    else if (!strcmp(argv[i],"-splitpbedge") ) 
      SplitEdgeWith2Boundary=1;
    else if (!strcmp(argv[i],"-nosplitpbedge") ) 
      SplitEdgeWith2Boundary=0;
    else if (!strcmp(argv[i],"-NbJacobi")  && ++i<argc) 
      nbjacoby = atoi(argv[i]);
    else if (!strcmp(argv[i],"-CutOff")  && ++i<argc) 
      CutOff =  atof(argv[i]), AbsError = 0;
    else if (!strcmp(argv[i],"-NbSmooth")  && ++i<argc) 
      NbSmooth = atoi(argv[i]);
    else if (!strcmp(argv[i],"-omega")  && ++i<argc) 
      omega =  atof(argv[i]);
    else 
      { 
	cout << " Error in the arguments " << i << "  of " << bamgversion <<endl;
	cout << "      argument [" << i << "]= `" << ((i< argc ) ?  argv[i] : " The Lastest ") 
	     <<"`" << endl;
	cout << " Usage :" << endl;
	cout << "  Mesh INPUT: The 2 arguments are exclusives." << endl;
	cout << "" << endl;
	cout << "     -g  filename    Set the geometry for mesh generation. " << endl;
	cout << "                     DB mesh file." << endl;
	cout << "     -b  filename    Set the background mesh for mesh adaption " << endl;
	cout << "                     (require {\tt -M} or {-Mbb} arguments). " << endl;
	cout << "                     - am_fmt file if filename match  *.am_fmt " << endl;
	cout << "                     - ambda file if filename match  *.amdba " << endl;
	cout << "                     - nopo file if filename match  *.nopo " << endl;
	cout << "                     - am file if filename match  *.am " << endl;
	cout << "                     - ftq file if filename match  *.ftq " << endl;
	cout << "                     - msh file if filename match  *.msh " << endl;
	cout << "                     - otherwise the file is a BD mesh file " << endl;
	cout << "                     Remark: the geometry is the background geometry." << endl;
	cout << "                     DB mesh file." << endl;
	cout << "     -r  filename    Set the  mesh for modification mesh  " << endl;
	cout << "                     no reconstruction " << endl;
	cout << "                     same as in the case of -b argument" << endl;
	cout << "" << endl;
	cout << "     -thetamax (double)   change the angular limit for a corner in degre " << endl;
	cout << "                       the angle is defined from 2 normals of 2 concecutives edges " << endl;
	cout << "                       if no geomtry cf reading an  am_fmt, amdba, nopo .. file " << endl;


	cout << "" << endl;
	cout << "  METRIC definition or mesh size definition, one of the" << endl;
	cout << "  2 next arguments is need in case of  mesh adaption." << endl;
	cout << "" << endl;
	cout << "     -M filename     Set  the metric which is  defined on the background mesh " << endl;
	cout << "                     or on the geometry. Metric file." << endl;
	cout << "     -Mbb filename   Set the solution  defined on the background mesh for" << endl;
	cout << "                     metric construction, the solutions was FE P1 defined on " << endl;
	cout << "                     the background mesh. bb file. " << endl;
	cout << "     -MBB filename   same with -Mbb but with BB file " << endl;
	cout << "" << endl;
	cout << "     -errg (double)  Set the level of error on geometry (0.1) by default" << endl;
	cout << "     -coef (double)  Set the value of mutiplicative " << endl;
        cout << "                     coef on the mesh size (1 by default)." << endl;
	cout << "     -power (double) Set the power parameter of hessien to construct " << endl;
	cout << "                     the metric  (1 by default)" << endl;
	cout << "     -maxsubdiv  (double) Change the metric  such that  the maximal subdivision  " << endl;
	cout << "                     of a background's edge is bound by the " <<endl;
	cout << "                     given number (always limited by 10) " << endl;
	cout << "     -ratio (double) Set the ratio for a smoothing of the metric." << endl;
	cout << "                     If ratio is  0 then no smoothing" << endl;
	cout << "	              and if ratio  is in  [1.1,10] then the smoothing " << endl;
	cout << "                     change the metrix such that the greatest geometrical" << endl;
	cout << "                     progression (speed of mesh size variation) " << endl;
	cout << "                     in a mesh is bounded  by ratio (by default no smoothing)." << endl;
	cout << "     -hmin (double)  Set the value of the minimal edge size." << endl;
	cout << "     -hminaniso (double)  Set the value of the minimal edge size and save aniso." << endl;
	cout << "     -hmax (double)  Set the value of the maximal edge size." << endl;
	cout << "     -NbSmooth (int) Number of Smoothing iteration " << endl;
	cout << "                     (3 by default if the metric is set otherwise 0)." << endl;
	cout << "     -omega (double)  relaxation parameter for Smoothing " << endl;
	cout << "     -splitpbedge     split in 2 all internal edges with 2 vertex on boundary" << endl;
	cout << "     -nosplitpbedge   d'ont cut internal edges with 2 vertex on boundary (default)" << endl;
	cout << "" << endl;
	cout << "" << endl;
	cout << "        the next arguments are used with the -Mbb argument" << endl;
	cout << "" << endl;
	cout << "     -KeepBackVertices  Try to Keep old vertices (default) " << endl;
	cout << "     -noKeepBackVertices  all vertices are create from scratch  " << endl;
	cout << "     -err   (double)    Set the level of error (default 0.01)"    << endl;
	cout << "     -iso               The constructed metric must be isotropic" << endl;
	cout << "     -aniso             The constructed metric can be anisotropic " << endl;
	cout << "     -anisomax (double) Set  maximum ratio  of anisotropy " << endl;
	cout << "                            1 =>  isotropic (1e6 by default) " << endl;
	cout << "     -RelError          Construct metric with relative  error " << endl;
	cout << "     -AbsError          Construct metric with with abs error " << endl;
	cout << "     -CutOff (double)   Set the limit of the relative error  (1.e-5)" << endl;
	cout << "     -NbJacobi (int)    Set the number of Jacobi iteration for smoothing " << endl;
	cout << "                        the construction of metric  (1 by default)." << endl;
	cout << "     -NoRescaling       Don't rescaling of all solution between [0..1] " 
	     << "before metric computation " << endl; 
	cout << "     -Rescaling       Do rescaling of all solution between [0..1] " 
	     << "before metric computation ((default case)" << endl; 
	cout << "     -H (int)           choices for computing the hessein (test)" << endl;
	cout << "" << endl;
	cout << "" << endl;
	cout << "  Definition of some internal variable and limitation." << endl;
	cout << "" << endl;
	cout << "     -v   (int)       Set the level of impression (verbosity) " << endl;
	cout << "                      between 0 and 10 (1 by default)." << endl;
	cout << "     -nbv (int)       Set the maximal of  vertices (" << nbvx  << " by default)." << endl;
	cout << "" << endl;
	cout << "" << endl;
	cout << "  To interpoled a solution form background mesh to generated " << endl;
	cout << "  mesh (in case of adpatation)" << endl;
	cout << "" << endl;
	cout << "    -rbb filename    Read  solution  file defined on the background " << endl;
	cout << "                      mesh for interpolation on created mesh." << endl;
	cout << "                      bb file. (by default the  -Mbb filename)" << endl;
	cout << "     -wbb filename    Write the file of interpolation of the solutions" << endl;
	cout << "                      read with  -rbb argument. " << endl;
	cout << "                      bb file." << endl;
	cout << "     -rBB filename    same -rbb but with BB file ";
	cout << "     -wBB filename    same -wbb but with BB file ";
	cout << "" << endl;
	cout << "  Output mesh file for adpation or generation." << endl;
	cout << "    Remark: one of output mesh file is require " << endl;
	cout << endl;
	cout << "     -o       filename Create a DB Mesh file." << endl;
	cout << "     -oamdba  filename Create a amdba file." << endl;
	cout << "     -oftq    filename Create a ftq file." << endl;
	cout << "     -omsh    filename Create a msh file (freefem3 file)." << endl;
	cout << "     -oam_fmt filename Create a am_fmt file. " << endl;
	cout << "     -oam     filename Create a am file. " << endl;
	cout << "     -onopo   filename Create a nopo file. " << endl;
	cout << "     -oM filename      Create a metric file. " << endl;

	cout << endl;
	cout << endl;
	cout << "     -thetaquad (double)  minimal angle of a quadrangle " << endl;
	cout << "     -2q                  split triangles in 3 quad and quad in 4 quad  " << endl;
	cout << "     -2                   split triangles in 4 trai and quad in 4 quad  " << endl;

	cout << endl;
      
	cout << " Remark: if no argument is given when the arguments are read in "<<datafile