ex4.php

一个用来实现偏微分方程中网格的计算库

      
      <B><FONT COLOR="#5F9EA0">std</FONT></B>::string order = <B><FONT COLOR="#BC8F8F">&quot;SECOND&quot;</FONT></B>; 
      <B><FONT COLOR="#A020F0">if</FONT></B> ( command_line.search(2, <B><FONT COLOR="#BC8F8F">&quot;-Order&quot;</FONT></B>, <B><FONT COLOR="#BC8F8F">&quot;-o&quot;</FONT></B>) )
        order = command_line.next(order);
  
      <B><FONT COLOR="#5F9EA0">std</FONT></B>::string family = <B><FONT COLOR="#BC8F8F">&quot;LAGRANGE&quot;</FONT></B>; 
      <B><FONT COLOR="#A020F0">if</FONT></B> ( command_line.search(2, <B><FONT COLOR="#BC8F8F">&quot;-FEFamily&quot;</FONT></B>, <B><FONT COLOR="#BC8F8F">&quot;-f&quot;</FONT></B>) )
        family = command_line.next(family);
      
      <B><FONT COLOR="#A020F0">if</FONT></B> ((family == <B><FONT COLOR="#BC8F8F">&quot;MONOMIAL&quot;</FONT></B>) || (family == <B><FONT COLOR="#BC8F8F">&quot;XYZ&quot;</FONT></B>))
        {
  	<B><FONT COLOR="#5F9EA0">std</FONT></B>::cerr &lt;&lt; <B><FONT COLOR="#BC8F8F">&quot;ex4 currently requires a C^0 (or higher) FE basis.&quot;</FONT></B> &lt;&lt; std::endl;
  	error();
        }
        
      Mesh mesh (dim);
      
  
      <B><FONT COLOR="#A020F0">if</FONT></B> ((family == <B><FONT COLOR="#BC8F8F">&quot;LAGRANGE&quot;</FONT></B>) &amp;&amp; (order == <B><FONT COLOR="#BC8F8F">&quot;FIRST&quot;</FONT></B>))
        {
  	<B><FONT COLOR="#5F9EA0">MeshTools</FONT></B>::Generation::build_cube (mesh,
  					   ps, ps, ps,
  					   -1., 1.,
  					   -1., 1.,
  					   -1., 1.,
  					   (dim==1)    ? EDGE2 : 
  					   ((dim == 2) ? QUAD4 : HEX8));
        }
      
      <B><FONT COLOR="#A020F0">else</FONT></B>
        {
  	<B><FONT COLOR="#5F9EA0">MeshTools</FONT></B>::Generation::build_cube (mesh,
  					   ps, ps, ps,
  					   -1., 1.,
  					   -1., 1.,
  					   -1., 1.,
  					   (dim==1)    ? EDGE3 : 
  					   ((dim == 2) ? QUAD9 : HEX27));
        }
  
      
      mesh.print_info();
      
      
      EquationSystems equation_systems (mesh);
      
      {
        LinearImplicitSystem&amp; system =
  	equation_systems.add_system&lt;LinearImplicitSystem&gt; (<B><FONT COLOR="#BC8F8F">&quot;Poisson&quot;</FONT></B>);
  
        
        system.add_variable(<B><FONT COLOR="#BC8F8F">&quot;u&quot;</FONT></B>,
  			  <B><FONT COLOR="#5F9EA0">Utility</FONT></B>::string_to_enum&lt;Order&gt;   (order),
  			  <B><FONT COLOR="#5F9EA0">Utility</FONT></B>::string_to_enum&lt;FEFamily&gt;(family));
  
        system.attach_assemble_function (assemble_poisson);
        
        equation_systems.init();
  
        equation_systems.print_info();
      }
  
      equation_systems.get_system(<B><FONT COLOR="#BC8F8F">&quot;Poisson&quot;</FONT></B>).solve();
  
      <B><FONT COLOR="#A020F0">if</FONT></B>(dim == 1)
      {        
        GnuPlotIO plot(mesh,<B><FONT COLOR="#BC8F8F">&quot;Example 4, 1D&quot;</FONT></B>,GnuPlotIO::GRID_ON);
        plot.write_equation_systems(<B><FONT COLOR="#BC8F8F">&quot;out_1&quot;</FONT></B>,equation_systems);
      }
      <B><FONT COLOR="#A020F0">else</FONT></B>
      {
        GMVIO (mesh).write_equation_systems ((dim == 3) ? 
          <B><FONT COLOR="#BC8F8F">&quot;out_3.gmv&quot;</FONT></B> : <B><FONT COLOR="#BC8F8F">&quot;out_2.gmv&quot;</FONT></B>,equation_systems);
      }
    }
    
    <B><FONT COLOR="#A020F0">return</FONT></B> libMesh::close ();
  }
  
  
  
  
  <B><FONT COLOR="#228B22">void</FONT></B> assemble_poisson(EquationSystems&amp; es,
                        <B><FONT COLOR="#228B22">const</FONT></B> std::string&amp; system_name)
  {
    assert (system_name == <B><FONT COLOR="#BC8F8F">&quot;Poisson&quot;</FONT></B>);
  
  
    PerfLog perf_log (<B><FONT COLOR="#BC8F8F">&quot;Matrix Assembly&quot;</FONT></B>);
    
    <B><FONT COLOR="#228B22">const</FONT></B> Mesh&amp; mesh = es.get_mesh();
  
    <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> dim = mesh.mesh_dimension();
  
    LinearImplicitSystem&amp; system = es.get_system&lt;LinearImplicitSystem&gt;(<B><FONT COLOR="#BC8F8F">&quot;Poisson&quot;</FONT></B>);
    
    <B><FONT COLOR="#228B22">const</FONT></B> DofMap&amp; dof_map = system.get_dof_map();
  
    FEType fe_type = dof_map.variable_type(0);
  
    AutoPtr&lt;FEBase&gt; fe (FEBase::build(dim, fe_type));
    
    QGauss qrule (dim, FIFTH);
  
    fe-&gt;attach_quadrature_rule (&amp;qrule);
  
    AutoPtr&lt;FEBase&gt; fe_face (FEBase::build(dim, fe_type));
  	      
    QGauss qface(dim-1, FIFTH);
    
    fe_face-&gt;attach_quadrature_rule (&amp;qface);
  
    <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;Real&gt;&amp; JxW = fe-&gt;get_JxW();
  
    <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;Point&gt;&amp; q_point = fe-&gt;get_xyz();
  
    <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;std::vector&lt;Real&gt; &gt;&amp; phi = fe-&gt;get_phi();
  
    <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;std::vector&lt;RealGradient&gt; &gt;&amp; dphi = fe-&gt;get_dphi();
  
    DenseMatrix&lt;Number&gt; Ke;
    DenseVector&lt;Number&gt; Fe;
  
    <B><FONT COLOR="#5F9EA0">std</FONT></B>::vector&lt;<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>&gt; dof_indices;
  
    <B><FONT COLOR="#5F9EA0">MeshBase</FONT></B>::const_element_iterator       el     = mesh.local_elements_begin();
    <B><FONT COLOR="#228B22">const</FONT></B> MeshBase::const_element_iterator end_el = mesh.local_elements_end();
  
    <B><FONT COLOR="#A020F0">for</FONT></B> ( ; el != end_el; ++el)
      {
        perf_log.start_event(<B><FONT COLOR="#BC8F8F">&quot;elem init&quot;</FONT></B>);      
  
        <B><FONT COLOR="#228B22">const</FONT></B> Elem* elem = *el;
  
        dof_map.dof_indices (elem, dof_indices);
  
        fe-&gt;reinit (elem);
  
        Ke.resize (dof_indices.size(),
  		 dof_indices.size());
  
        Fe.resize (dof_indices.size());
  
        perf_log.stop_event(<B><FONT COLOR="#BC8F8F">&quot;elem init&quot;</FONT></B>);      
  
        perf_log.start_event (<B><FONT COLOR="#BC8F8F">&quot;Ke&quot;</FONT></B>);
  
        <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> qp=0; qp&lt;qrule.n_points(); qp++)
  	<B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i&lt;phi.size(); i++)
  	  <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> j=0; j&lt;phi.size(); j++)
  	    Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
  	    
  
        perf_log.stop_event (<B><FONT COLOR="#BC8F8F">&quot;Ke&quot;</FONT></B>);
  
        perf_log.start_event (<B><FONT COLOR="#BC8F8F">&quot;Fe&quot;</FONT></B>);
        
        <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> qp=0; qp&lt;qrule.n_points(); qp++)
  	{
  	  <B><FONT COLOR="#228B22">const</FONT></B> Real x = q_point[qp](0);
  	  <B><FONT COLOR="#228B22">const</FONT></B> Real y = q_point[qp](1);
  	  <B><FONT COLOR="#228B22">const</FONT></B> Real z = q_point[qp](2);
  	  <B><FONT COLOR="#228B22">const</FONT></B> Real eps = 1.e-3;
  
  	  <B><FONT COLOR="#228B22">const</FONT></B> Real uxx = (exact_solution(x-eps,y,z) +
  			    exact_solution(x+eps,y,z) +
  			    -2.*exact_solution(x,y,z))/eps/eps;
  	      
  	  <B><FONT COLOR="#228B22">const</FONT></B> Real uyy = (exact_solution(x,y-eps,z) +
  			    exact_solution(x,y+eps,z) +
  			    -2.*exact_solution(x,y,z))/eps/eps;
  	  
  	  <B><FONT COLOR="#228B22">const</FONT></B> Real uzz = (exact_solution(x,y,z-eps) +
  			    exact_solution(x,y,z+eps) +
  			    -2.*exact_solution(x,y,z))/eps/eps;
  
            Real fxy;
            <B><FONT COLOR="#A020F0">if</FONT></B>(dim==1)
            {
              <B><FONT COLOR="#228B22">const</FONT></B> Real pi = libMesh::pi;
              fxy = (0.25*pi*pi)*sin(.5*pi*x);
            }
            <B><FONT COLOR="#A020F0">else</FONT></B>
            {
  	    fxy = - (uxx + uyy + ((dim==2) ? 0. : uzz));
            } 
  
  	  <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i&lt;phi.size(); i++)
  	    Fe(i) += JxW[qp]*fxy*phi[i][qp];	  
  	}
        
        perf_log.stop_event (<B><FONT COLOR="#BC8F8F">&quot;Fe&quot;</FONT></B>);
  
        {
  	
  	perf_log.start_event (<B><FONT COLOR="#BC8F8F">&quot;BCs&quot;</FONT></B>);
  
  	<B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> side=0; side&lt;elem-&gt;n_sides(); side++)
  	  <B><FONT COLOR="#A020F0">if</FONT></B> (elem-&gt;neighbor(side) == NULL)
  	    {
              
  	      <B><FONT COLOR="#228B22">const</FONT></B> Real penalty = 1.e10;
  
                <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;std::vector&lt;Real&gt; &gt;&amp;  phi_face = fe_face-&gt;get_phi();
  
                <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;Real&gt;&amp; JxW_face = fe_face-&gt;get_JxW();
  
                <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;Point &gt;&amp; qface_point = fe_face-&gt;get_xyz();
  
                fe_face-&gt;reinit(elem, side);
  
                <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> qp=0; qp&lt;qface.n_points(); qp++)
                {
                  <B><FONT COLOR="#228B22">const</FONT></B> Real xf = qface_point[qp](0);
                  <B><FONT COLOR="#228B22">const</FONT></B> Real yf = qface_point[qp](1);
                  <B><FONT COLOR="#228B22">const</FONT></B> Real zf = qface_point[qp](2);
  
  
                  <B><FONT COLOR="#228B22">const</FONT></B> Real value = exact_solution(xf, yf, zf);
  
                  <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i&lt;phi_face.size(); i++)
                    <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> j=0; j&lt;phi_face.size(); j++)
                      Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
  
                  <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i&lt;phi_face.size(); i++)
                    Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp];
                } 
              }
              
  	
  	perf_log.stop_event (<B><FONT COLOR="#BC8F8F">&quot;BCs&quot;</FONT></B>);
        } 
        
  
        perf_log.start_event (<B><FONT COLOR="#BC8F8F">&quot;matrix insertion&quot;</FONT></B>);
        
        system.matrix-&gt;add_matrix (Ke, dof_indices);
        system.rhs-&gt;add_vector    (Fe, dof_indices);
  
        perf_log.stop_event (<B><FONT COLOR="#BC8F8F">&quot;matrix insertion&quot;</FONT></B>);
      }
  
  }
</pre> 
<a name="output"></a> 
<br><br><br> <h1> The console output of the program: </h1> 
<pre>
***************************************************************
* Running Example  ./ex4-devel
***************************************************************
 
Running ./ex4-devel -d 1 -n 20

 Mesh Information:
  mesh_dimension()=1
  spatial_dimension()=3
  n_nodes()=41
  n_elem()=20
   n_local_elem()=20
   n_active_elem()=20
  n_subdomains()=1
  n_processors()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System "Poisson"
    Type "LinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE" 
    Approximation Orders="SECOND" 
    n_dofs()=41
    n_local_dofs()=41
    n_constrained_dofs()=0
    n_vectors()=1


-------------------------------------------------------
| Time:           Wed Jun  6 12:18:07 2007             |
| OS:             Linux                                |
| HostName:       orville                              |
| OS Release:     2.6.21-1.3194.fc7PAE                 |
| OS Version:     #1 SMP Wed May 23 22:27:31 EDT 2007  |
| Machine:        i686                                 |
| Username:       peterson                             |
-------------------------------------------------------
 ------------------------------------------------------------------------------
| Matrix Assembly Performance: Alive time=0.001533, Active time=0.000777       |
 ------------------------------------------------------------------------------
| Event                         nCalls    Total       Avg         Percent of   |
|                                         Time        Time        Active Time  |
|------------------------------------------------------------------------------|
|                                                                              |
| BCs                           20        0.0001      0.000005    13.26        |
| Fe                            20        0.0002      0.000010    26.25        |
| Ke                            20        0.0000      0.000002    5.53         |
| elem init                     20        0.0002      0.000011    28.44        |
| matrix insertion              20        0.0002      0.000010    26.51        |
 ------------------------------------------------------------------------------
| Totals:                       100       0.0008                  100.00       |
 ------------------------------------------------------------------------------

Running ./ex4-devel -d 2 -n 15

 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=961
  n_elem()=225
   n_local_elem()=225
   n_active_elem()=225
  n_subdomains()=1
  n_processors()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System "Poisson"
    Type "LinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE" 
    Approximation Orders="SECOND" 
    n_dofs()=961
    n_local_dofs()=961
    n_constrained_dofs()=0
    n_vectors()=1


-------------------------------------------------------
| Time:           Wed Jun  6 12:18:07 2007             |
| OS:             Linux                                |
| HostName:       orville                              |
| OS Release:     2.6.21-1.3194.fc7PAE                 |
| OS Version:     #1 SMP Wed May 23 22:27:31 EDT 2007  |
| Machine:        i686                                 |
| Username:       peterson                             |
-------------------------------------------------------
 ------------------------------------------------------------------------------
| Matrix Assembly Performance: Alive time=0.02087, Active time=0.018112        |
 ------------------------------------------------------------------------------
| Event                         nCalls    Total       Avg         Percent of   |
|                                         Time        Time        Active Time  |
|------------------------------------------------------------------------------|
|                                                                              |
| BCs                           225       0.0042      0.000019    23.18        |
| Fe                            225       0.0046      0.000021    25.52        |
| Ke                            225       0.0038      0.000017    21.18        |
| elem init                     225       0.0037      0.000016    20.30        |
| matrix insertion              225       0.0018      0.000008    9.82         |
 ------------------------------------------------------------------------------
| Totals:                       1125      0.0181                  100.00       |
 ------------------------------------------------------------------------------

Running ./ex4-devel -d 3 -n 6

 Mesh Information:
  mesh_dimension()=3
  spatial_dimension()=3
  n_nodes()=2197
  n_elem()=216
   n_local_elem()=216
   n_active_elem()=216
  n_subdomains()=1
  n_processors()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System "Poisson"
    Type "LinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE" 
    Approximation Orders="SECOND" 
    n_dofs()=2197
    n_local_dofs()=2197
    n_constrained_dofs()=0
    n_vectors()=1


-------------------------------------------------------
| Time:           Wed Jun  6 12:18:07 2007             |
| OS:             Linux                                |
| HostName:       orville                              |
| OS Release:     2.6.21-1.3194.fc7PAE                 |
| OS Version:     #1 SMP Wed May 23 22:27:31 EDT 2007  |
| Machine:        i686                                 |
| Username:       peterson                             |
-------------------------------------------------------
 ------------------------------------------------------------------------------
| Matrix Assembly Performance: Alive time=0.240115, Active time=0.237185       |
 ------------------------------------------------------------------------------
| Event                         nCalls    Total       Avg         Percent of   |
|                                         Time        Time        Active Time  |
|------------------------------------------------------------------------------|
|                                                                              |
| BCs                           216       0.1129      0.000523    47.59        |
| Fe                            216       0.0128      0.000059    5.41         |
| Ke                            216       0.0812      0.000376    34.22        |
| elem init                     216       0.0169      0.000078    7.12         |
| matrix insertion              216       0.0134      0.000062    5.66         |
 ------------------------------------------------------------------------------
| Totals:                       1080      0.2372                  100.00       |
 ------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example  ./ex4-devel
***************************************************************
</pre>
</div>
<?php make_footer() ?>
</body>
</html>
<?php if (0) { ?>
\#Local Variables:
\#mode: html
\#End:
<?php } ?>