medit.cpp

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

// ORIG-DATE:     Aout 2008
// -*- Mode : c++ -*-
//
// SUMMARY  : liaison medit freefem++ : popen  
// USAGE    : LGPL      
// ORG      : LJLL Universite Pierre et Marie Curie, Paris,  FRANCE 
// AUTHOR   : Jacques Morice
// E-MAIL   : [e-mail]
//

/* 
 This file is part of 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

 Thank to the ARN ()  FF2A3 grant
 ref:ANR-07-CIS7-002-01 
 */
#include "mode_open.hpp"
#include <iostream>
#include <cfloat>
#include <cmath>
#include <complex>
using namespace std;
#include "error.hpp"
#include "AFunction.hpp"
using namespace std;  
#include "rgraph.hpp"
#include "RNM.hpp"
#include "fem.hpp"


#include "FESpacen.hpp" 
#include "FESpace.hpp" 

#include "MatriceCreuse_tpl.hpp"
#include "MeshPoint.hpp"
#include "Operator.hpp" 
#include "lex.hpp"

#include "lgfem.hpp"
#include "lgmesh3.hpp"
#include "libmesh5.h"
#include "lgsolver.hpp"
#include "problem.hpp"
//#include "LayerMesh.hpp"
//#include "TransfoMesh_v2.hpp"
//#include "GQuadTree.hpp"

#include <set>
#include <vector>
#include <fstream>

//#include "Operator.hpp" 
//#include "array_resize.hpp"
//#include "lex.hpp"

//using Fem2D::Mesh;
//using Fem2D::MeshPoint;

//extern bool NoWait; 
//const char *medit_char= "meditff"; *
#define WrdSiz 4

#ifdef WIN32
string stringffmedit= "ffmedit.exe";
//string stringemptymedit= "ffmedit.exe";
#else
string stringffmedit= "ffmedit";
#endif


const char *medit_popen="-popen";// 1";  // depend de l endroit ou se trouve medit
const char *medit_bin="-filebin";
const char *medit_addsol="-addsol";

const char *medit_debug="-d";

static bool TheWait=false;
bool  NoWait=false;

using namespace std;
using namespace Fem2D;


//*************************
//
//   creation point sol
//
//*************************
// datasolMesh2

class datasolMesh2_Op : public E_F0mps 
{
public:
  typedef long  Result;
  Expression eTh;
  Expression filename;
  
  struct Expression2 {
    long what; // 1 scalar, 2 vector, 3 symtensor
    long nbfloat; // 1 scalar, 2 vector (3D), 3 symtensor(3D)
    Expression e[3];
    Expression2() {e[0]=0; e[1]=0; e[2]=0;  what=0; nbfloat=0;};
    Expression &operator[](int i){return e[i];}
    double eval(int i,Stack stack) const  { 
    if (e[i]) {
      return GetAny< double >( (*e[i])(stack) );
    }
    else 
      return 0;
    }
  };
  vector<Expression2> l;
  static const int n_name_param = 1;  
  static basicAC_F0::name_and_type name_param[];
  Expression nargs[n_name_param];
  long    arg(int i,Stack stack,long a ) const{ return nargs[i] ? GetAny< long >( (*nargs[i])(stack) ): a;}
 
public:
  datasolMesh2_Op(const basicAC_F0 &  args) : l( args.size()-2 )
  {
    int nbofsol;
    int ddim=2;
    int stsize=3;
    //cout << "construction data medit solution avec datasolMesh2_Op" << args << endl;
    //cout << "taille de args" << args.size() << endl;
    args.SetNameParam(n_name_param,name_param,nargs); 

    if (BCastTo<string *>(args[0])) filename = CastTo<string *>(args[0]);
    if (BCastTo<pmesh>(args[1])) eTh= CastTo<pmesh>(args[1]);
   
  
    nbofsol = l.size();
    for (size_t i=2;i<args.size();i++){
      size_t jj=i-2;

      if ( BCastTo<double>( args[i] ))
	{
	  l[jj].what=1;
	  l[jj].nbfloat=1;
	  l[jj][0]=to<double>( args[i] );
	  
	}
      else if ( args[i].left()==atype<E_Array>() )
	{
	  const E_Array * a0  = dynamic_cast<const E_Array *>( args[i].LeftValue() );
	  //cout << "taille" << a0->size() << endl;
	  //if (a0->size() != ddim || a0->size() != stsize) 
	  //  CompileError("savesol in 2D: vector solution is 2 composant, vector solution is 3 composant");
	  
	  if( a0->size() == ddim){
	    // vector solution
	    l[jj].what=2;
	    l[jj].nbfloat=ddim;
	    for(int j=0; j<ddim; j++){
	      l[jj][j] = to<double>( (*a0)[j]);
	    }
	  }
	  else if( a0->size() == stsize){
	    // symmetric tensor solution
	    l[jj].what=3;
	    l[jj].nbfloat=stsize;
	    for(int j=0; j<stsize; j++){
	      l[jj][j] = to<double>( (*a0)[j]);
	    }
	  }
	  
	}
      else {
	cout << " arg " << i << " " << args[i].left() << endl;
	CompileError("savesol in 2D: Sorry no way to save this kind of data");
      }
      
    }
  }
  static ArrayOfaType  typeargs() { return  ArrayOfaType( atype<string *>(), atype<pmesh>(), true); }// all type
  static  E_F0 * f(const basicAC_F0 & args) { return new datasolMesh2_Op(args);} 
  AnyType operator()(Stack stack)  const ;
};

basicAC_F0::name_and_type datasolMesh2_Op::name_param[]= {
  {  "order",&typeid(long)}
};

AnyType datasolMesh2_Op::operator()(Stack stack)  const 
{ 
  MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
  Mesh * pTh= GetAny<Mesh *>((*eTh)(stack));
  string * ffname= GetAny<string *>( (*filename)(stack) );
  ffassert(pTh);
  Mesh &Th=*pTh;

  int nt = Th.nt;
  int nv = Th.nv;

  int nbtype=l.size();
  int nbsol;
  int solnbfloat;
  int TypTab[l.size()];
 
  int resultorder= arg(0, stack, 1);
  long longdefault;

  int ver = GmfFloat, outm;
  // determination de TypTab
  solnbfloat=0;
  for (size_t i=0;i<l.size();i++){
    TypTab[i]=l[i].what;
    solnbfloat=solnbfloat+l[i].nbfloat;
  }
  float *OutSolTab = new float[solnbfloat];

  // determination de OutSolTab

  if ( !(outm = GmfOpenMesh(ffname->c_str(),GmfWrite,ver,2)) ) {
    cerr <<"  -- Mesh3::Save  UNABLE TO OPEN  :"<< ffname << endl;
    exit(1);
  }
 
  if(resultorder==0){
    // ordre 0
    nbsol = nt;  
    
    KN<double> valsol(solnbfloat*nbsol);
 
    MeshPoint *mp3(MeshPointStack(stack)); 
    R2 Cdg_hat = R2(1./3.,1./3.);  

    for (int it=0;it<nt;it++){
      int h=0;
      const Mesh::Triangle  & K(Th.t(it));
      mp3->set( Th, K(Cdg_hat), Cdg_hat, K, K.lab);
	  
      for(size_t i=0;i<l.size();i++){
	for(size_t j=0;j<l[i].nbfloat;j++){
	  valsol[it*solnbfloat+h] = l[i].eval(j,stack);
	  h=h+1;
	}
      } 
      assert(solnbfloat==h);
    }


    GmfSetKwd(outm,GmfSolAtTriangles, nbsol, nbtype, TypTab);
    for (int k=0; k<nbsol; k++){
      for (int i=0; i<solnbfloat ;i++){
	OutSolTab[i] =  valsol(k*solnbfloat+i);
      }
      GmfSetLin(outm, GmfSolAtTriangles, OutSolTab); 
    }
  }
  if(resultorder==1){
    // ordre 1
    nbsol = nv;  

    KN<double> valsol(solnbfloat*nbsol);
    valsol=0.;
    KN<int> takemesh(nbsol);
    MeshPoint *mp3(MeshPointStack(stack)); 
    //R2 Cdg_hat = R2(1./3.,1./3.);
    takemesh=0;
    for (int it=0;it<nt;it++){
      for(int iv=0;iv<3;iv++){
	int i=Th(it,iv);

	//if(takemesh[i]==0){
	mp3->setP(&Th,it,iv);
	int h=0;
	
	for(size_t ii=0;ii<l.size();ii++){
	  for(size_t j=0;j<l[ii].nbfloat;j++){
	    //cout << "ii=" << ii << " j=" << j<< endl;
	    valsol[i*solnbfloat+h] = valsol[i*solnbfloat+h] + l[ii].eval(j,stack);
	    h=h+1;
	  }
	} 
	assert(solnbfloat==h);
	takemesh[i] = takemesh[i]+1;
	//}
      }
    }
    for(int i=0; i<nv; i++){
      for(int h=0; h<solnbfloat; h++){
	valsol[i*solnbfloat+h] = valsol[i*solnbfloat+h]/takemesh[i]; 
      }
    }

    GmfSetKwd(outm,GmfSolAtVertices, nbsol, nbtype, TypTab);
    for (int k=0; k<nbsol; k++){
      for (int i=0; i<solnbfloat ;i++){
	OutSolTab[i] =  valsol(k*solnbfloat+i);
      }
      GmfSetLin(outm, GmfSolAtVertices, OutSolTab);
    }
  }
  GmfCloseMesh(outm);

  delete [] OutSolTab;
  return longdefault;
}
  
// datasolMesh3
template<class v_fes>
class datasolMesh3_Op : public E_F0mps 
{
public:
  typedef long  Result;
  Expression eTh;
  Expression filename;
  
  struct Expression2 {
    long what; // 1 scalar, 2 vector, 3 symtensor
    long nbfloat; // 1 scalar, 3 vector (3D), 6 symtensor(3D)
    Expression e[6];
    Expression2() {e[0]=0; e[1]=0; e[2]=0; e[3]=0; e[4]=0; e[5]=0; what=0; nbfloat=0;};
    Expression &operator[](int i){return e[i];}
    double eval(int i,Stack stack) const  { 
    if (e[i]) {
      return GetAny< double >( (*e[i])(stack) );
    }
    else 
      return 0;
    }
    
  };
  vector<Expression2> l; 
  static const int n_name_param =1;  
  static basicAC_F0::name_and_type name_param[];
  Expression nargs[n_name_param];
  long    arg(int i,Stack stack,long a ) const{ return nargs[i] ? GetAny< long >( (*nargs[i])(stack) ): a;}
 
public:
  datasolMesh3_Op(const basicAC_F0 &  args) : l( args.size()-2 )
  {
    int nbofsol;
    int ddim=3;
    int stsize=6;
    //cout << "construction data medit solution avec datasolMesh3_Op" << args << endl;
    args.SetNameParam(n_name_param,name_param,nargs); 

    //if (BCastTo<string *>(args[0])){
      filename = CastTo<string *>(args[0]);
      // }
      //if (BCastTo<pmesh3>(args[1])){
      eTh= CastTo<pmesh3>(args[1]);
      // }
    nbofsol = l.size();
    for (size_t i=2;i<args.size();i++){
      size_t jj=i-2;

      if ( BCastTo<double>(args[i]))
	{
	  l[jj].what=1;
	  l[jj].nbfloat=1;
	  l[jj][0]=to<double>( args[i] );
	  
	}
      else if ( args[i].left()==atype<E_Array>() )
	{
	  const E_Array * a0  = dynamic_cast<const E_Array *>( args[i].LeftValue() );
	  //cout << "taille" << a0->size() << endl;
	  //if (a0->size() != ddim || a0->size() != stsize) 
	  //  CompileError("savesol in 2D: vector solution is 2 composant, vector solution is 3 composant");

	  if( a0->size() == ddim){
	    // vector solution
	    l[jj].what=2;
	    l[jj].nbfloat=ddim;
	    for(int j=0; j<ddim; j++){
	      l[jj][j] = to<double>( (*a0)[j]);
	    }
	  }
	  else if( a0->size() == stsize){
	    // symmetric tensor solution
	    l[jj].what=3;
	    l[jj].nbfloat=stsize;
	    for(int j=0; j<stsize; j++){
	      l[jj][j] = to<double>( (*a0)[j]);
	    }
	  }
	  
	}
      else {
	CompileError("savesol in 3D: Sorry no way to save this kind of data");
      }

    }
  }
  static ArrayOfaType  typeargs() { return  ArrayOfaType( atype<string *>(), atype<pmesh3>(), true); }// all type
  static  E_F0 * f(const basicAC_F0 & args) { return new datasolMesh3_Op(args);} 
  AnyType operator()(Stack stack)  const ;
};
template<class v_fes>
basicAC_F0::name_and_type datasolMesh3_Op<v_fes>::name_param[]= {
  { "order",&typeid(long)}
};
template<class v_fes>
AnyType datasolMesh3_Op<v_fes>::operator()(Stack stack)  const 
{ 
  MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
  Mesh3 * pTh= GetAny<Mesh3 *>((*eTh)(stack));
  string * ffname= GetAny<string *>( (*filename)(stack) );
  ffassert(pTh);
  Mesh3 &Th=*pTh;

  int nt = Th.nt;
  int nv = Th.nv;
  int nbe = Th.nbe;
  int nbtype=l.size();
  int nbsol;
  int solnbfloat;
  int TypTab[l.size()];
 
  int resultorder= arg(0, stack, 1);
  long longdefault;

  int ver = GmfFloat, outm;
  // determination de TypTab
  solnbfloat=0;
  for (size_t i=0;i<l.size();i++){
    TypTab[i]=l[i].what;
    solnbfloat=solnbfloat+l[i].nbfloat;
  }
  float *OutSolTab = new float[solnbfloat];

  // determination de OutSolTab

  if ( !(outm = GmfOpenMesh(ffname->c_str(),GmfWrite,ver,3)) ) {
    cerr <<"  -- Mesh3::Save  UNABLE TO OPEN  :"<< filename << endl;
    exit(1);
  }
 
  if(resultorder==0){
    // Tetrahedra
    // ordre 0
    nbsol = nt;      
    KN<double> valsol(solnbfloat*nbsol);
 
    MeshPoint *mp3(MeshPointStack(stack)); 
    R3 Cdg_hat = R3(1./4.,1./4.,1./4.);  

    for (int it=0;it<nt;it++){
      int h=0;
      const Tet & K(Th.elements[it]);
      mp3->set( Th, K(Cdg_hat), Cdg_hat, K, K.lab);
	  
      for(size_t i=0;i<l.size();i++){
	for(size_t j=0;j<l[i].nbfloat;j++){
	  valsol[it*solnbfloat+h] = l[i].eval(j,stack);
	  h=h+1;
	}
      } 
      assert(solnbfloat==h);
    }
    
    GmfSetKwd(outm,GmfSolAtTetrahedra, nbsol, nbtype, TypTab);
    for (int k=0; k<nbsol; k++){
      for (int i=0; i<solnbfloat ;i++){
	OutSolTab[i] =  valsol(k*solnbfloat+i);
      }
      GmfSetLin(outm, GmfSolAtTetrahedra, OutSolTab); 
    }
    
    /*
    // Rajout des Triangles
    nbsol = nbe;      
    
    KN<double> valsol2(solnbfloat*nbsol);
    //cout << "value of triangles in the border nbe="<< nbe << endl;
    MeshPoint *mp32(MeshPointStack(stack)); 
     
    for (int it=0;it<nbe;it++){
      int h=0;
      const Triangle3 & K(Th.be(it));
      R xg,yg,zg;
      int iv[3];
      for(int jj=0; jj <3; jj++){
	iv[jj] = Th.operator()(K[jj]); 
      }

      xg = Th.vertices[iv[0]].x + Th.vertices[iv[1]].x + Th.vertices[iv[2]].x;
      yg = Th.vertices[iv[0]].y + Th.vertices[iv[1]].y + Th.vertices[iv[2]].y;
      zg = Th.vertices[iv[0]].z + Th.vertices[iv[1]].z + Th.vertices[iv[2]].z;

      xg=xg/3.;
      yg=yg/3.;
      zg=zg/3.;

      mp32->set(xg,yg,zg);
	  
      for(size_t i=0;i<l.size();i++){
	for(size_t j=0;j<l[i].nbfloat;j++){
	  valsol2[it*solnbfloat+h] = l[i].eval(j,stack);
	  h=h+1;
	}
      } 
      assert(solnbfloat==h);
    }

    GmfSetKwd(outm,GmfSolAtTriangles, nbsol, nbtype, TypTab);
    for (int k=0; k<nbsol; k++){
      for (int i=0; i<solnbfloat ;i++){
	OutSolTab[i] =  valsol2(k*solnbfloat+i);
      }
      //cout <<"GmfSetLin(outm, GmfSolAtTriangles, OutSolTab);"<< endl;
      GmfSetLin(outm, GmfSolAtTriangles, OutSolTab); 
    }

    */
  }
  if(resultorder==1){
    // ordre 1
    nbsol = nv;  

    KN<double> valsol(solnbfloat*nbsol);
    KN<int> takemesh(nbsol);
    MeshPoint *mp3(MeshPointStack(stack)); 
    R3 Cdg_hat = R3(1./4.,1./4.,1./4.);
    takemesh=0;
    for (int it=0;it<nt;it++){
      for(int iv=0;iv<4;iv++){
	int i=Th(it,iv);

	if(takemesh[i]==0){
	    mp3->setP(&Th,it,iv);
	    int h=0;