drawing.cpp

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

// -*- Mode : c++ -*-
//
// SUMMARY  :      
// USAGE    :        
// ORG      : 
// AUTHOR   : Frederic Hecht
// 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
 */

#include <cmath>
#include "error.hpp"
#include <cstdio>
#include <iostream>
#include <iomanip>
#include <fstream>
#include "rgraph.hpp"


#include "RNM.hpp"
#include "fem.hpp"
#include "FESpacen.hpp" 
#include "FESpace.hpp" 

namespace Fem2D {

void NewSetColorTable(int nb,float *colors=0,int nbcolors=0,bool hsv=true)
{
  if(colors && nbcolors)
    SetColorTable1(nb,hsv,nbcolors,colors);
  else
    SetColorTable(nb);
}
int dichotomie(const RN_ &viso,R v) 
{
  int i=0,j=viso.N(),k;
  if  (v <viso[0] || v >viso[j-1]) 
    return -1;  
  while (i<j-1)    
   if ( viso[k=(i+j)/2]> v) j=k;
   else i=k;
  return i;
}

void plot(long i)
{
  char buf[24];
  snprintf(buf,24,"%ld",i);
  plotstring(buf);
}
void plot(int i)
{
  char buf[24];
  snprintf(buf,24,"%d",i);
  plotstring(buf);
}
void plot(double i)
{
  char buf[24];
  snprintf(buf,24,"%g",i);
  plotstring(buf);
}

void DrawIsoT(const R2 Pt[3],const R ff[3],const RN_ & Viso);
void DrawIsoTfill(const R2 Pt[3],const R ff[3],const RN_ & Viso,double rapz);

void FillRect(float x0,float y0, float x1, float y1)
 {
     float r[8];
     r[0]=x0;r[1]=y0;
     r[2]=x1;r[3]=y0;
     r[4]=x1;r[5]=y1;
     r[6]=x0;r[7]=y1;
     fillpoly(4,r);
 }

void PlotValue(const RN_ & Viso,int  k0,const char * cmm)
{
   float xmin,xmax,ymin,ymax;
 //  cout << "PlotValue " << Viso << endl;
   // int ix,iy;
  // GetSizeScreen(ix,iy);   
   getcadre(xmin,xmax,ymin,ymax);
   float dx=(xmax-xmin);
   float dy=(ymax-ymin);
   //  10 points  
  // int kk = Max(30,iy/10);
   float h=GetHeigthFont();
   float x0=xmin+dx*0.8;
   float y= ymin+dy*0.95;
  // cout << x0 << " " << y << " " << h <<  endl;
   y -= k0* h*1.4;
   couleur(0);
   FillRect(x0-h*0.5,y-h*(1.4*Viso.N()+0.3),x0+h*9,y+h*1.5);
   couleur(1);
   rmoveto(x0+1.2*h,y);
   plotstring(cmm);
   y -=  h*1.4;
   for (int i=0;i<Viso.N();i++)
    {
       couleur(i+4);
       FillRect(x0,y,x0+h,y+h);
       rmoveto(x0+1.2*h,y);
       y -=  h*1.4;
       plot(Viso[i]);                 
    }
}

void DrawCommentaire(const char * cm,float x,float y) 
{
   float xmin,xmax,ymin,ymax;
   getcadre(xmin,xmax,ymin,ymax);
   float dx=(xmax-xmin);
   float dy=(ymax-ymin);
   
   rmoveto(xmin+dx*x,ymin+dy*y);
   plotstring(cm);   
}

void Mesh::DrawBoundary() const
{  
    for(int i=0;i<neb;i++)
       bedges[i].Draw();
}
void Mesh::InitDraw() const 
   {
     R2 Pmin,Pmax;
     BoundingBox(Pmin,Pmax);
     R2 O((Pmin+Pmax)/2);
     R r =  (Max(Pmax.x-Pmin.x,Pmax.y-Pmin.y)*0.55);
     showgraphic();
     cadreortho((float)O.x,(float)(O.y+r*0.05),(float) r);
    }
inline  R Square(R x){return x*x;}

void  SetDefaultIsoValue(const RN_& U,RN_ & Viso)
 {
   R umn = U.min();
   R umx = U.max();
   int N = Viso.N();
   R d = (umx-umn)/(N+1);
   R x = umn+d/2;
   for (int i = 0;i < N;i++)
     {Viso[i]=x;x +=d; }
   NewSetColorTable(N+4) ;  
 }
void  SetDefaultIsoValue(const RN_& u,const RN_& v,RN_ & Viso)
 {
   R s = (Square(u[0])+Square(v[0]));
   R umn =s ;
   R umx = s;    
   int n(v.N());
    for (int i=1;i<n;i++)
     {   s = (Square(u[i])+Square(v[i]));
        umn = Min(umn,s),umx=Max(umn,s);
      }  
   umx=sqrt(umx);
   umn=sqrt(umn);
   
   int N = Viso.N();
   R d = (umx-umn)/(N-1);
   R x = umn;
   for (int i = 0;i < N;i++)
     {Viso[i]=x;x +=d; }
   Viso[N-1]= umx; 
   NewSetColorTable(N+4) ;  
 }
 
 void Mesh::Draw(int init,bool fill) const
 {
   
   if (init>0)  InitDraw();
   
   SetColorTable(16) ;  
   for (int i=0;i<nt;i++)
    {
      if (fill) 
       {
         triangles[i].Fill(2+Abs(triangles[i].lab));         
         couleur(triangles[i].lab?1:0);
         triangles[i].Draw();
       }
      else 
       {
       couleur(1+Abs(triangles[i].lab));
       triangles[i].Draw();}
    }
     

   couleur(1);
   penthickness(2); 
   for (int i=0;i<NbMortars;i++)
    {     
      mortars[i].Draw();
    }
   penthickness(3);
   double l= 1e100;
   for (int i=0;i<neb;i++)
     {
       couleur(1+Abs( bedges[i].lab));
       bedges[i].Draw();
       l=Min(l,bedges[i].length());
     }
   for (int i=0;i<nv;i++)
     { 
       R2 P=vertices[i], N=vertices[i].Ne();
       if (Norme2_2(N) ) 
         { MoveTo(P); LineTo(P+N*l/5);  }
     }
   penthickness(1); 
   if(init%2)
   for (int i=0;i<nv;i++)
    {
      couleur(1+Abs(vertices[i].lab));
      MoveTo(vertices[i]);
      plot(i);
    }
   couleur(1);
//   rattente(0);
 }

void FElement::Draw(const RN_& U,const RN_ & Viso,int composante) const
{   
  int nsb = nbsubdivision();
  int nsb2 = NbOfSubTriangle(nsb);
  int nbdf=NbDoF();
  RN fk(nbdf);
  for (int i=0;i<nbdf;i++) // get the local value
    fk[i] = U[operator()(i)];
  RNMK fb(nbdf,N,3); //  the value for basic fonction
  RN ff(3);
  R2 Pt,P[3],A(T[0]),B(T[1]),C(T[2]);
  bool whatd[last_operatortype];
  initwhatd(whatd,0);
  
  for (int k=0;k<nsb2;k++)
   {   
   //  ff=0.0;
     for (int j=0;j<3;j++)
      { 
     //   cout << nbdf << endl;
        Pt=SubTriangle(nsb,k,j); //  current point 
        BF(whatd,Pt,fb);
        ff[j] = (fb('.',composante,0),fk);
        P[j] = A*(1-Pt.x-Pt.y)+B*Pt.x + C*Pt.y;
       }
     DrawIsoT(P,ff,Viso);  
   } 
} 
    
void FElement::Drawfill(const RN_& U,const RN_ & Viso,int composante, double rapz) const
{   
  int nsb = nbsubdivision();
  int nsb2 =  NbOfSubTriangle(nsb);
  int nbdf=NbDoF();
  RN fk(nbdf);
  for (int i=0;i<nbdf;i++) // get the local value
    fk[i] = U[operator()(i)];
  RNMK fb(nbdf,N,3); //  the value for basic fonction
  RN ff(3);
  R2 Pt,P[3],A(T[0]),B(T[1]),C(T[2]);
  bool whatd[last_operatortype];
  initwhatd(whatd,0);
  
  for (int k=0;k<nsb2;k++)
   {   
   //  ff=0.0;
     for (int j=0;j<3;j++)
      { 
     //   cout << nbdf << endl;
        Pt=SubTriangle(nsb,k,j); //  current point 
        BF(whatd,Pt,fb);
        ff[j] = (fb('.',composante,0),fk);
        P[j] = A*(1-Pt.x-Pt.y)+B*Pt.x + C*Pt.y;
       }
     DrawIsoTfill(P,ff,Viso,rapz);  
   } 
} 

R2 FElement::MinMax(const RN_& U,const RN_& V,int i0,int i1) const
{
  R2 minmax(1e100,-1e100);
  int nsb = nbsubdivision();
  int nsb2 =  NbOfSubTriangle(nsb);
  int nbdf=NbDoF();
  RN fk(nbdf);
  RN gk(nbdf);
  RNMK fb(nbdf,N,3); //  the value for basic fonction
  RN f0(3);
  RN f1(3);
  bool whatd[last_operatortype];
  initwhatd(whatd,0);  
  R2 Pt,A(T[0]),B(T[1]),C(T[2]);
  for (int i=0;i<nbdf;i++) // get the local value
    fk[i] = U[operator()(i)];
  for (int i=0;i<nbdf;i++) // get the local value
    gk[i] = V[operator()(i)];
  for (int k=0;k<nsb2;k++)
   {   
   //  ff=0.0;
     for (int j=0;j<3;j++)
      { 
     //   cout << nbdf << endl;
        R2 Pt=SubTriangle(nsb,k,j); //  current point 
        BF(whatd,Pt,fb);
        f0[j] = (fb('.',i0,0),fk);
        f1[j] = (fb('.',i1,0),gk);
        R2 uv(f0[j],f1[j]);
        R uv2 = (uv,uv);
        minmax.x=Min(minmax.x,uv2);
        minmax.y=Max(minmax.y,uv2);
        
       }
    }
   return minmax;
}

R2 FElement::MinMax(const RN_& U,int i0) const
{
  R2 minmax(1e100,-1e100);
  int nsb = nbsubdivision();
  int nsb2 =  NbOfSubTriangle(nsb);
  int nbdf=NbDoF();
  RN fk(nbdf);
  RNMK fb(nbdf,N,3); //  the value for basic fonction
  RN f0(3);
  RN f1(3);
  R2 Pt,A(T[0]),B(T[1]),C(T[2]);
   bool whatd[last_operatortype];
  initwhatd(whatd,0);
 
  for (int i=0;i<nbdf;i++) // get the local value
    fk[i] = U[operator()(i)];
  for (int k=0;k<nsb2;k++)
   {   
   //  ff=0.0;
     for (int j=0;j<3;j++)
      { 
     //   cout << nbdf << endl;
        R2 Pt=SubTriangle(nsb,k,j); //  current point 
        BF(whatd,Pt,fb);
        f0[j] = (fb('.',i0,0),fk);
        minmax.x=Min(minmax.x,f0[j]);
        minmax.y=Max(minmax.y,f0[j]);        
       }
    }
   return minmax;
}
void FElement::Draw(const RN_& U,const RN_& V,const RN_ & Viso,R coef,int i0,int i1) const
{
  bool whatd[last_operatortype];
  initwhatd(whatd,0);
   
  float xmin,xmax,ymin,ymax;
  getcadre(xmin,xmax,ymin,ymax);
  float d= Max(ymax-ymin,xmax-xmin);
  R kk = d*0.005;
  R cc = d*0.05;
  int nsb = nbsubdivision();
  int nsb2 = NbOfSubTriangle(nsb);
  int nbdf=NbDoF();
  RN fk(nbdf);