-- FreeFem++ v4.6 (Thu Apr  2 15:47:38 CEST 2020 - git v4.6)
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : // 3D parallel plate capacitor problem
    2 : //
    3 : //  Discussion:
    4 : //
    5 : //    Solves the Laplace Equation on a parallel plate capacitor (dirichlet) boundary.
    6 : //    FE-Mesh and space functions are written to text files in order to be plot with
    7 : //    Matlab / Octave.
    8 : //
    9 : //    This program requires the ffmatlib.idp include file, as well as access to
   10 : //    the msh3 and medit libraries.
   11 : //
   12 : //  Licensing:
   13 : //
   14 : //    Copyright (C) 2018 Chloros2 <chloros2@gmx.de>
   15 : //
   16 : //    This program is free software: you can redistribute it and/or modify it
   17 : //    under the terms of the GNU General Public License as published by
   18 : //    the Free Software Foundation, either version 3 of the License, or
   19 : //    (at your option) any later version.
   20 : //
   21 : //    This program is distributed in the hopeC that it will be useful, but
   22 : //    WITHOUT ANY WARRANTY; without even the implied warranty of
   23 : //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   24 : //    GNU General Public License for more details.
   25 : //
   26 : //    You should have received a copy of the GNU General Public License
   27 : //    along with this program.  If not, see
   28 : //    <https://www.gnu.org/licenses/>.
   29 : //
   30 : //  Modified:
   31 : //
   32 : //    2018-05-19
   33 : //
   34 : // Author: 
   35 : //
   36 : //   Chloros2 <chloros2@gmx.de>
   37 : //
   38 : include "ffmatlib.idp"/*ffmatlib.idp Export Macros (Plot Solutions for Matlab and Octave)
    2 :  *
    3 :  * Author: Chloros2 <chloros2@gmx.de>
    4 :  * Created: 2018-12-21
    5 :  *
    6 :  *   This file is part of the ffmatlib which is hosted at
    7 :  *   https://github.com/samplemaker/freefem_matlab_octave_plot
    8 :  *
    9 :  * Copyright (C) 2018 Chloros2 <chloros2@gmx.de>
   10 :  *
   11 :  * This program is free software: you can redistribute it and/or modify it
   12 :  * under the terms of the GNU General Public License as published by
   13 :  * the Free Software Foundation, either version 3 of the License, or
   14 :  * (at your option) any later version.
   15 :  *
   16 :  * This program is distributed in the hope that it will be useful, but
   17 :  * WITHOUT ANY WARRANTY; without even the implied warranty of
   18 :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   19 :  * GNU General Public License for more details.
   20 :  *
   21 :  * You should have received a copy of the GNU General Public License
   22 :  * along with this program.  If not, see
   23 :  * <https://www.gnu.org/licenses/>.
   24 :  *
   25 :  */
   26 : 
   27 : //Save finite element connectivity (j-th degree of freedom for i-th element)
   28 : macro ffSaveVh(Th, Vh, filename){
   29 #   {
   30 #   ofstream file(filename);
   31 #   file.precision(10);
   32 #   for (int i=0; i<Th.nt; i++){
   33 #     for (int j=0; j<Vh.ndofK; j++){
   34 #       file << Vh(i,j) << "\n";
   35 #     }
   36 #   }
   37 #   }
   38 # }  )  //EOM
   39 : 
   40 : //Save FE Space data
   41 : macro ffSaveData(u, filename){
   42 #   {
   43 #   ofstream file(filename);
   44 #   file.precision(10);
   45 #   for (int j=0; j<u[].n; j++){
   46 #     file << u[][j] << endl;
   47 #   }
   48 #   }
   49 # }  )  //EOM
   50 : 
   51 : macro ffSaveData2(u1, u2, filename){
   52 #   {
   53 #   ofstream file(filename);
   54 #   int datalen=u1[].n;
   55 #   if (u2[].n!=datalen){
   56 #     cout << "error: arguments must have same size" << endl;
   57 #     exit(1);
   58 #   }
   59 #   file.precision(10);
   60 #   for (int j=0; j<datalen; j++){
   61 #     file << u1[][j] << " " << u2[][j] << endl;
   62 #   }
   63 #   }
   64 # }  )  //EOM
   65 : 
   66 : macro ffSaveData3(u1, u2, u3, filename){
   67 #   {
   68 #   ofstream file(filename);
   69 #   int datalen=u1[].n;
   70 #   if ((u2[].n!=datalen) | (u3[].n!=datalen)){
   71 #     cout << "error: arguments must have same size" << endl;
   72 #     exit(1);
   73 #   }
   74 #   file.precision(10);
   75 #   for (int j=0; j<datalen; j++){
   76 #     file << u1[][j] << " " << u2[][j] << " " << u3[][j] << endl;
   77 #   }
   78 #   }
   79 # }  )  //EOM
   80 : 
   81 : macro ffSaveData4(u1, u2, u3, u4, filename){
   82 #   {
   83 #   ofstream file(filename);
   84 #   int datalen=u1[].n;
   85 #   if ((u2[].n!=datalen) | (u3[].n!=datalen) | (u4[].n!=datalen)){
   86 #     cout << "error: arguments must have same size" << endl;
   87 #     exit(1);
   88 #   }
   89 #   file.precision(10);
   90 #   for (int j=0; j<datalen; j++){
   91 #     file << u1[][j] << " " << u2[][j] << " " << u3[][j] << " " << u4[][j] << endl;
   92 #   }
   93 #   }
   94 # }  )  //EOM
   95 : 
   39 : 
   40 : cout << endl;
   41 : cout << "capacitor_3d:" << endl;
   42 : cout << "  FreeFem++ version." << endl;
   43 : cout << "  Model a capacitor." << endl;
   44 : cout << "  Using the ffmatlib interface, export the data f
  ... : or graphic" << endl;
   45 : cout << "  processing by MATLAB/Octave." << endl;
   46 : 
   47 : load "msh3"
   48 : load "medit"
   49 : 
   50 : int C20=98,C21=99; //2D mesh
   51 : int CK=30,CA=31;   //anode+cathode
   52 : int C3CT01=23, C3CT12=24, C3CT23=25, C3CT34=26; //3D contact labels
   53 : int C3BOT=40, C3TOP=41, C3SIDEO=42, C3SIDEI=43; //3D labels - cube top and bottom and side labels
   54 : 
   55 : //enclosure
   56 : real w0=0.75, h0=0.75;
   57 : 
   58 : border C01(t=0,w0){ x=t; y=0; label=C20;};
   59 : border C02(t=w0,0){ x=t; y=h0; label=C20;};
   60 : border C03(t=0,h0){ x=w0; y=t; label=C20;};
   61 : border C04(t=h0,0){ x=0; y=t; label=C20;};
   62 : 
   63 : real w1=0.3, h1=0.3;
   64 : real posX1=0.5*w0-0.5*w1;
   65 : real posY1=0.5*h0-0.5*h1;
   66 : 
   67 : border C11(t=0,w1){ x=posX1+t; y=posY1; label=C21;};
   68 : border C12(t=w1,0){ x=posX1+t; y=posY1+h1; label=C21;};
   69 : border C13(t=0,h1){ x=posX1+w1; y=posY1+t; label=C21;};
   70 : border C14(t=h1,0){ x=posX1; y=posY1+t; label=C21;};
   71 : 
   72 : //layer stack
   73 : real[int] z(6)
 *** Warning  The identifier z hide a Global identifier  
;
   74 : z=[0.0,0.17,0.2,0.3,0.33,0.5];
   75 : 
   76 : int n=30;
   77 : int m=12;
   78 : real l=22.0/(z(5)-z(0));
   79 : 
   80 : mesh Tha2=buildmesh(C01(n)+C03(n)+C02(n)+C04(n)
   81 :                    +C11(m)+C12(m)+C13(m)+C14(m));
   82 : 
   83 : //bottom block
   84 : int[int] a0up=[0,CK,          //Deckel Innen
   85 :                1,C3CT01],     //Deckel Aussenring
   86 :          a0down=[0,C3BOT,     //Boden Innen
   87 :                  1,C3BOT],    //Boden Aussenring
   88 :          a0mid=[C20,C3SIDEO,  //Mantel Aussen
   89 :                 C21,C3SIDEI], //Mantel Innen
   90 :          a0tet=[0,0];
   91 : 
   92 : int[int] a1up=[0,CK,
   93 :                1,C3CT12],
   94 :          a1down=[0,CK,
   95 :                  1,C3CT01],
   96 :          a1mid=[C20,C3SIDEO,
   97 :                 C21,CK],
   98 :          a1tet=[0,0];
   99 : 
  100 : int[int] a2up=[0,CA,
  101 :                1,C3CT23],
  102 :          a2down=[0,CK,
  103 :                  1,C3CT12],
  104 :          a2mid=[C20,C3SIDEO,
  105 :                 C21,C3SIDEI],
  106 :          a2tet=[0,0];
  107 : 
  108 : int[int] a3up=[0,CA,
  109 :                1,C3CT34],
  110 :          a3down=[0,CA,
  111 :                  1,C3CT23],
  112 :          a3mid=[C20,C3SIDEO,
  113 :                 C21,CA],
  114 :          a3tet=[0,0];
  115 : 
  116 : int[int] a4up=[0,55,
  117 :                1,C3TOP],
  118 :          a4down=[0,CA,
  119 :                  1,C3CT34],
  120 :          a4mid=[C20,C3SIDEO,
  121 :                 C21,C3SIDEI],
  122 :          a4tet=[0,0];
  123 : 
  124 : mesh3 Tha03d=buildlayers(Tha2,l*(z[1]-z[0]),
  125 :                          zbound=[z[0],z[1]],
  126 :                          reftet=a0tet,
  127 :                          labelmid=a0mid,
  128 :                          labelup=a0up,
  129 :                          labeldown=a0down);
  130 : 
  131 : mesh3 Tha13d=buildlayers(Tha2,l*(z[2]-z[1]),
  132 :                          zbound=[z[1],z[2]],
  133 :                          reftet=a1tet,
  134 :                          labelmid=a1mid,
  135 :                          labelup=a1up,
  136 :                          labeldown=a1down);
  137 : 
  138 : mesh3 Tha23d=buildlayers(Tha2,l*(z[3]-z[2]),
  139 :                          zbound=[z[2],z[3]],
  140 :                          reftet=a2tet,
  141 :                          labelmid=a2mid,
  142 :                          labelup=a2up,
  143 :                          labeldown=a2down);
  144 : 
  145 : mesh3 Tha33d=buildlayers(Tha2,l*(z[4]-z[3]),
  146 :                          zbound=[z[3],z[4]],
  147 :                          reftet=a3tet,
  148 :                          labelmid=a3mid,
  149 :                          labelup=a3up,
  150 :                          labeldown=a3down);
  151 : 
  152 : mesh3 Tha43d=buildlayers(Tha2,l*(z[5]-z[4]),
  153 :                          zbound=[z[4],z[5]],
  154 :                          reftet=a4tet,
  155 :                          labelmid=a4mid,
  156 :                          labelup=a4up,
  157 :                          labeldown=a4down);
  158 : 
  159 : mesh3 Thn3d=Tha03d+Tha13d+Tha23d+Tha33d+Tha43d;
  160 : 
  161 : fespace Vh(Thn3d,P1);
  162 : Vh u,v,h;
  163 : 
  164 : //Terminal K on Ground
  165 : //Arbitrary Anode Voltage
  166 : real u0=2.0;
  167 : 
  168 : macro Grad(u) [dx(u),dy(u),dz(u)]  )  // EOM
  169 : 
  170 : problem Laplace(u,v,solver=CG) =
  171 :                int3d(Thn3d)( Grad(u)   [dx(u),dy(u),dz(u)] '*Grad(v)   [dx(v),dy(v),dz(v)]  )
  172 :                +on(CK,u=0)
  173 :                +on(CA,u=u0);
  174 : 
  175 : Laplace;
  176 : 
  177 : Vh Ex, Ey, Ez;
  178 : Ex = -dx(u);
  179 : Ey = -dy(u);
  180 : Ez = -dz(u);
  181 : 
  182 : real epsilon0 = 8.854187e-12; // [As/Vm]
  183 : real energy = 0.5*epsilon0*int3d(Thn3d,qforder=3)( (Ex)^2 + (Ey)^2 + (Ez)^2 );
  184 : 
  185 : cout.precision(3);
  186 : 
  187 : cout << endl;
  188 : cout << "A Parallel Plate Capacitor Problem:" << endl;
  189 : cout << "Capacitance:	" << "C = " << 2*energy/(u0)^2 << "[F]" << endl;
  190 : cout << "epsilon*A/d:	" << "C = " << epsilon0*w1*h1/(z[3]-z[2]) << "[F]" << endl;
  191 : cout << endl;
  192 : cout << "NbBoundaryElements:	" <<  Thn3d.nbe << endl;
  193 : cout << "NbTriangles:		" <<  Thn3d.nt << endl;
  194 : cout << "NbVertices:		" <<  Thn3d.nv << endl;
  195 : cout << "nDoF:			" << Vh.ndof << endl;
  196 : cout << endl;
  197 : //
  198 : //  Write solution data to files for processing by MATLAB.
  199 : //
  200 : savemesh (Thn3d, "capacitor_3d.mesh" );
  201 : ffSaveVh (Thn3d, Vh, "capacitor_3d_vh.txt" )
   29 @   
   30 @    
   31 @   
   32 @       
   33 @         
   34 @           
   35 @     
   36 @   
   37 @   
   38 @  {
   29 @   {
   30 @   ofstream file( "capacitor_3d_vh.txt");
   31 @   file.precision(10);
   32 @   for (int i=0; i<Thn3d.nt; i++){
   33 @     for (int j=0; j< Vh.ndofK; j++){
   34 @       file <<  Vh(i,j) << "\n";
   35 @     }
   36 @   }
   37 @   }
   38 @ } ;
  202 : ffSaveData ( u, "capacitor_3d_pot.txt" )
   42 @   
   43 @    
   44 @   
   45 @       
   46 @         
   47 @   
   48 @   
   49 @  {
   42 @   {
   43 @   ofstream file( "capacitor_3d_pot.txt");
   44 @   file.precision(10);
   45 @   for (int j=0; j< u[].n; j++){
   46 @     file <<  u[][j] << endl;
   47 @   }
   48 @   }
   49 @ } ;
  203 : ffSaveData3 ( Ex, Ey, Ez, "capacitor_3d_vec.txt" )
   67 @   
   68 @    
   69 @    
   70 @      
   71 @         
   72 @     
   73 @   
   74 @   
   75 @       
   76 @                 
   77 @   
   78 @   
   79 @  {
   67 @   {
   68 @   ofstream file( "capacitor_3d_vec.txt");
   69 @   int datalen= Ex[].n;
   70 @   if (( Ey[].n!=datalen) | ( Ez[].n!=datalen)){
   71 @     cout << "error: arguments must have same size" << endl;
   72 @     exit(1);
   73 @   }
   74 @   file.precision(10);
   75 @   for (int j=0; j<datalen; j++){
   76 @     file <<  Ex[][j] << " " <<  Ey[][j] << " " <<  Ez[][j] << endl;
   77 @   }
   78 @   }
   79 @ } ;
  204 : 
  205 : medit("U",Thn3d,u);
  206 : //
  207 : //  Terminate.
  208 : //
  209 : cout << "\n";
  210 : cout << "capacitor_3d:\n";
  211 : cout << "  Normal end of execution.\n";
  212 : 
  213 : exit ( 0 );
  214 : 
  215 :  sizestack + 1024 =4784  ( 3760 )


capacitor_3d:
  FreeFem++ version.
  Model a capacitor.
  Using the ffmatlib interface, export the data for graphic
  processing by MATLAB/Octave.
  --  mesh:  Nb of Triangles =   2018, Nb of Vertices 1070
  -- FESpace: Nb of Nodes 22470 Nb of DoF 22470
 GC:  converge in  164 g=7.05392e-12 rho= 1.8066 gamma= 0.766907
  -- Solve : 
          min 1.24523e-33  max 2
  -- Medit,  Release 3.0a (Nov. 30, 2007)
     Copyright (c) LJLL, 1999-2007.
     compiled:  Thu Apr 2 15:52:17 CEST 2020 (with ff++ 4.6).

medit with binary version of popen : mesh(es) and solution(s) 
mesh_name= U
 Loading data file(s)
End of mesh
.sol: Dimension 3 (mesh)3 (lecture)1 
  Input seconds:     0.02

 medit1() 

 Building scene(s)
  Creating scene 1
   Loading default options
  Scene seconds:     0.05

 Rendering scene(s)

 Total running seconds:  0.12
 Thank you for using Medit.

A Parallel Plate Capacitor Problem:
Capacitance:	C = 1.54e-11[F]
epsilon*A/d:	C = 7.97e-12[F]

NbBoundaryElements:	18828
NbTriangles:		121080
NbVertices:		22470
nDoF:			22470

version de medit ffmedit -popen -filebin -addsol 1  U

capacitor_3d:
  Normal end of execution.
  current line = 213
exit(0)
 err code 0 ,  mpirank 0
 CodeAlloc : nb ptr  4519,  size :512088 mpirank: 0
Ok: Normal End