msh3.cpp

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

// ORIG-DATE: Novembre 2008
// -*- Mode : c++ -*-
//
// SUMMARY  :  
// 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 
 */
#ifndef WITH_NO_INIT
/*
#include <fstream>
#include <iostream>
#include <cfloat>
#include <cmath>
#include <cstring>
#include <complex>
//#include <set>
#include<stdlib.h>
//#include <vector>
#include <map>
using namespace std;
#include "error.hpp"
#include "AFunction.hpp"
#include "ufunction.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 "Mesh2dn.hpp"
#include "Mesh3dn.hpp"
#include "Operator.hpp" 
#include "lex.hpp"
#include "libmesh5.h"
#include "lgfem.hpp"
#include "lgmesh3.hpp"
#include "lgsolver.hpp"
#include "problem.hpp"
*/
#include "ff++.hpp"
#endif

//  TransfoMesh_v2.cpp
using namespace std;
// LayerMesh.cpp
// buildlayer.cpp
// rajout global

#include <set>
#include <vector>
#include "msh3.hpp"

using namespace  Fem2D;

void TestSameVertexMesh3( const Mesh3 & Th3, const double & hseuil, const R3 & Psup, const R3 &Pinf, int & nv_t, int *Numero_Som){
  
  Vertex3  *v=new Vertex3[Th3.nv];
  nv_t=0;
  
  EF23::GTree<Vertex3> *gtree = new EF23::GTree<Vertex3>(v,Pinf,Psup,0);
  
  // creation of octree
  for (int ii=0;ii<Th3.nv;ii++){
    const R3 r3vi( Th3.vertices[ii].x, Th3.vertices[ii].y, Th3.vertices[ii].z );
    const Vertex3 &vi(r3vi);
    
    Vertex3 * pvi=gtree->ToClose(vi,hseuil);
    
    if(!pvi){
      v[nv_t].x = vi.x;
      v[nv_t].y = vi.y;
      v[nv_t].z = vi.z;
      v[nv_t].lab = Th3.vertices[ii].lab; // lab mis a zero par default
      Numero_Som[ii] = nv_t; 
      gtree->Add( v[nv_t] );
      nv_t=nv_t+1;
    }
    else{
      Numero_Som[ii] = pvi-v;
    }
  }
  
  delete gtree;
  delete [] v;
}

void TestSameTetrahedraMesh3( const Mesh3 & Th3, const double & hseuil, const R3 & Psup, const R3 &Pinf, int & nt_t ){

  Vertex3 *vt=new Vertex3[Th3.nt];
  EF23::GTree<Vertex3> *gtree_t = new EF23::GTree<Vertex3>(vt,Pinf,Psup,0);
  
  nt_t=0;
  // creation of octree
  for (int ii=0;ii<Th3.nt;ii++){
    const Tet & K(Th3.elements[ii]);
    int iv[4];

    for(int jj=0; jj <4; jj++)
      iv[jj] = Th3.operator()(K[jj]);
    
    const double Cdg_x = ( Th3.vertices[iv[0]].x + Th3.vertices[iv[1]].x + Th3.vertices[iv[2]].x + Th3.vertices[iv[3]].x )/4.;
    const double Cdg_y = ( Th3.vertices[iv[0]].y + Th3.vertices[iv[1]].y + Th3.vertices[iv[2]].y + Th3.vertices[iv[3]].y )/4.;
    const double Cdg_z = ( Th3.vertices[iv[0]].z + Th3.vertices[iv[1]].z + Th3.vertices[iv[2]].z + Th3.vertices[iv[3]].z )/4.;

    const R3 r3vi( Cdg_x, Cdg_y, Cdg_z );
    const Vertex3 &vi(r3vi);
    
    Vertex3 * pvi=gtree_t->ToClose(vi,hseuil);
    
    if(!pvi){
      vt[nt_t].x = vi.x;
      vt[nt_t].y = vi.y;
      vt[nt_t].z = vi.z;
      vt[nt_t].lab = K.lab ; // lab mis a zero par default
      gtree_t->Add( vt[nt_t] );
      nt_t=nt_t+1;
    }
  }
  
  delete gtree_t;
  delete [] vt;
} 

void TestSameTetrahedraMesh3( const Mesh3 & Th3, const double & hseuil, const R3 & Psup, const R3 &Pinf, int *Elem_ok, int & nt_t ){

  Vertex3 *vt=new Vertex3[Th3.nt];
  EF23::GTree<Vertex3> *gtree_t = new EF23::GTree<Vertex3>(vt,Pinf,Psup,0);
  
  nt_t=0;
  // creation of octree
  for (int ii=0;ii<Th3.nt;ii++){
    if(Elem_ok[ii]!=1) continue;
    const Tet & K(Th3.elements[ii]);
    int iv[4];

    for(int jj=0; jj <4; jj++)
      iv[jj] = Th3.operator()(K[jj]);
    
    const double Cdg_x = ( Th3.vertices[iv[0]].x + Th3.vertices[iv[1]].x + Th3.vertices[iv[2]].x + Th3.vertices[iv[3]].x )/4.;
    const double Cdg_y = ( Th3.vertices[iv[0]].y + Th3.vertices[iv[1]].y + Th3.vertices[iv[2]].y + Th3.vertices[iv[3]].y )/4.;
    const double Cdg_z = ( Th3.vertices[iv[0]].z + Th3.vertices[iv[1]].z + Th3.vertices[iv[2]].z + Th3.vertices[iv[3]].z )/4.;

    const R3 r3vi( Cdg_x, Cdg_y, Cdg_z );
    const Vertex3 &vi(r3vi);
    
    Vertex3 * pvi=gtree_t->ToClose(vi,hseuil);
    
    if(!pvi){
      vt[nt_t].x = vi.x;
      vt[nt_t].y = vi.y;
      vt[nt_t].z = vi.z;
      vt[nt_t].lab = K.lab ; // lab mis a zero par default
      gtree_t->Add( vt[nt_t] );
      nt_t=nt_t+1;
    }
    else{
      Elem_ok[ii]=0;
    }
  }
  
  delete gtree_t;
  delete [] vt;
} 


void TestSameTriangleMesh3( const Mesh3 & Th3, const double & hseuil, const R3 & Psup, const R3 &Pinf, int & nbe_t){

  Vertex3  *vbe= new Vertex3[Th3.nbe];
  EF23::GTree<Vertex3> *gtree_be = new EF23::GTree<Vertex3>(vbe,Pinf,Psup,0);
  
  nbe_t=0;
  // creation of octree
  for (int ii=0;ii<Th3.nbe;ii++){
    const Triangle3 & K(Th3.be(ii));
    int iv[3];

    for(int jj=0; jj<3; jj++)
      iv[jj] = Th3.operator()(K[jj]);
    
    const double Cdg_x = ( Th3.vertices[iv[0]].x + Th3.vertices[iv[1]].x + Th3.vertices[iv[2]].x )/3.;
    const double Cdg_y = ( Th3.vertices[iv[0]].y + Th3.vertices[iv[1]].y + Th3.vertices[iv[2]].y )/3.;
    const double Cdg_z = ( Th3.vertices[iv[0]].z + Th3.vertices[iv[1]].z + Th3.vertices[iv[2]].z )/3.;

    const R3 r3vi( Cdg_x, Cdg_y, Cdg_z );
    const Vertex3 &vi(r3vi);
    
    Vertex3 * pvi=gtree_be->ToClose(vi,hseuil);
    
    if(!pvi){
      vbe[nbe_t].x = vi.x;
      vbe[nbe_t].y = vi.y;
      vbe[nbe_t].z = vi.z;
      vbe[nbe_t].lab = K.lab ; // lab mis a zero par default
      gtree_be->Add( vbe[nbe_t] );
      nbe_t=nbe_t+1;
    }
   
  }
  
  delete gtree_be;
  delete [] vbe;
} 

void TestSameTriangleMesh3( const Mesh3 & Th3, const double & hseuil, const R3 & Psup, const R3 &Pinf, int *Border_ok ,int & nbe_t ){

  Vertex3 *vbe=new Vertex3 [Th3.nbe];
  EF23::GTree<Vertex3> *gtree_be = new EF23::GTree<Vertex3>(vbe,Pinf,Psup,0);
  
  nbe_t=0;
  // creation of octree
  for (int ii=0;ii<Th3.nbe;ii++){
    if(Border_ok[ii]!=1) continue;
    const Triangle3 & K(Th3.be(ii));
    int iv[3];

    for(int jj=0; jj<3; jj++)
      iv[jj] = Th3.operator()(K[jj]);
    
    const double Cdg_x = ( Th3.vertices[iv[0]].x + Th3.vertices[iv[1]].x + Th3.vertices[iv[2]].x )/3.;
    const double Cdg_y = ( Th3.vertices[iv[0]].y + Th3.vertices[iv[1]].y + Th3.vertices[iv[2]].y )/3.;
    const double Cdg_z = ( Th3.vertices[iv[0]].z + Th3.vertices[iv[1]].z + Th3.vertices[iv[2]].z )/3.;

    const R3 r3vi( Cdg_x, Cdg_y, Cdg_z );
    const Vertex3 &vi(r3vi);
    
    Vertex3 * pvi=gtree_be->ToClose(vi,hseuil);
    
    if(!pvi){
      vbe[nbe_t].x = vi.x;
      vbe[nbe_t].y = vi.y;
      vbe[nbe_t].z = vi.z;
      vbe[nbe_t].lab = K.lab ; // lab mis a zero par default
      gtree_be->Add( vbe[nbe_t] );
      nbe_t=nbe_t+1;
    }
    else{      
      if(K.lab == vbe[pvi-vbe].lab )  Border_ok[ii] = 0;
    }
  }
  
  delete gtree_be;
  delete [] vbe;
} 

int TestElementMesh3( const Mesh3 & Th3 ) 
// Test si le maillage