ex7.php

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

<div class ="fragment">
<pre>
            f_system.set_frequencies_by_steps (frequency_in/n_frequencies,
        				       frequency_in,
        				       n_frequencies);
            
</pre>
</div>
<div class = "comment">
Use the parameters of the equation systems object to
tell the frequency system about the wave velocity and fluid
density.  The frequency system provides default values, but
these may be overridden, as shown here.
</div>

<div class ="fragment">
<pre>
            equation_systems.parameters.set&lt;Real&gt; ("wave speed") = 1.;
            equation_systems.parameters.set&lt;Real&gt; ("rho")        = 1.;
            
</pre>
</div>
<div class = "comment">
Initialize the data structures for the equation system.  <i>Always</i>
prior to this, the frequencies have to be communicated to the system.
</div>

<div class ="fragment">
<pre>
            equation_systems.init ();
            
</pre>
</div>
<div class = "comment">
Prints information about the system to the screen.
</div>

<div class ="fragment">
<pre>
            equation_systems.print_info ();
        
            for (unsigned int n=0; n &lt; n_frequencies; n++)
              {
</pre>
</div>
<div class = "comment">
Solve the system "Helmholtz" for the n-th frequency.  
Since we attached an assemble() function to the system,
the mass, damping and stiffness contributions will only
be assembled once.  Then, the system is solved for the
given frequencies.  Note that solve() may also solve 
the system only for specific frequencies.
</div>

<div class ="fragment">
<pre>
                f_system.solve (n,n);
        	
</pre>
</div>
<div class = "comment">
After solving the system, write the solution
to a GMV-formatted plot file, for every frequency.  
Now this is nice ;-) : we have the <i>identical</i> 
interface to the mesh write method as in the real-only 
case, but we output the real and imaginary 
part, and the magnitude, where the variable 
"p" is prepended with "r_", "i_", and "a_", 
respectively.
</div>

<div class ="fragment">
<pre>
                char buf[14];
        	sprintf (buf, "out%04d.gmv", n);
        	GMVIO(mesh).write_equation_systems (buf,
        					    equation_systems);
              }
            
</pre>
</div>
<div class = "comment">
Alternatively, the whole EquationSystems object can be
written to disk.  By default, the additional vectors are also
saved.
</div>

<div class ="fragment">
<pre>
            equation_systems.write ("eqn_sys.dat", libMeshEnums::WRITE);
          }
          
</pre>
</div>
<div class = "comment">
All done.  
</div>

<div class ="fragment">
<pre>
          return libMesh::close ();
        
        #endif 
        }
        
        
</pre>
</div>
<div class = "comment">
Here we define the matrix assembly routine for
the Helmholtz system.  This function will be
called to form the stiffness matrix and right-hand side.
</div>

<div class ="fragment">
<pre>
        void assemble_helmholtz(EquationSystems& es,
        			const std::string& system_name)
        {
        #ifdef USE_COMPLEX_NUMBERS
            
</pre>
</div>
<div class = "comment">
It is a good idea to make sure we are assembling
the proper system.
</div>

<div class ="fragment">
<pre>
          assert (system_name == "Helmholtz");
          
</pre>
</div>
<div class = "comment">
Get a constant reference to the mesh object.
</div>

<div class ="fragment">
<pre>
          const Mesh& mesh = es.get_mesh();
          
</pre>
</div>
<div class = "comment">
The dimension that we are in
</div>

<div class ="fragment">
<pre>
          const unsigned int dim = mesh.mesh_dimension();
          
</pre>
</div>
<div class = "comment">
Get a reference to our system, as before
</div>

<div class ="fragment">
<pre>
          FrequencySystem & f_system =
            es.get_system&lt;FrequencySystem&gt; (system_name);
          
</pre>
</div>
<div class = "comment">
A const reference to the DofMap object for this system.  The DofMap
object handles the index translation from node and element numbers
to degree of freedom numbers.
</div>

<div class ="fragment">
<pre>
          const DofMap& dof_map = f_system.get_dof_map();
          
</pre>
</div>
<div class = "comment">
Get a constant reference to the Finite Element type
for the first (and only) variable in the system.
</div>

<div class ="fragment">
<pre>
          const FEType& fe_type = dof_map.variable_type(0);
        
</pre>
</div>
<div class = "comment">
For the admittance boundary condition,
get the fluid density
</div>

<div class ="fragment">
<pre>
          const Real rho = es.parameters.get&lt;Real&gt;("rho");
          
</pre>
</div>
<div class = "comment">
In here, we will add the element matrices to the
<i>additional</i> matrices "stiffness_mass", "damping",
and the additional vector "rhs", not to the members 
"matrix" and "rhs".  Therefore, get writable
references to them
</div>

<div class ="fragment">
<pre>
          SparseMatrix&lt;Number&gt;&   stiffness      = f_system.get_matrix("stiffness");
          SparseMatrix&lt;Number&gt;&   damping        = f_system.get_matrix("damping");
          SparseMatrix&lt;Number&gt;&   mass           = f_system.get_matrix("mass");
          NumericVector&lt;Number&gt;&  freq_indep_rhs = f_system.get_vector("rhs");
          
</pre>
</div>
<div class = "comment">
Some solver packages (PETSc) are especially picky about
allocating sparsity structure and truly assigning values
to this structure.  Namely, matrix additions, as performed
later, exhibit acceptable performance only for identical
sparsity structures.  Therefore, explicitly zero the
values in the collective matrix, so that matrix additions
encounter identical sparsity structures.
</div>

<div class ="fragment">
<pre>
          SparseMatrix&lt;Number&gt;&  matrix           = *f_system.matrix;
          
</pre>
</div>
<div class = "comment">
------------------------------------------------------------------
Finite Element related stuff

<br><br>Build a Finite Element object of the specified type.  Since the
FEBase::build() member dynamically creates memory we will
store the object as an AutoPtr<FEBase>.  This can be thought
of as a pointer that will clean up after itself.
</div>

<div class ="fragment">
<pre>
          AutoPtr&lt;FEBase&gt; fe (FEBase::build(dim, fe_type));
          
</pre>
</div>
<div class = "comment">
A 5th order Gauss quadrature rule for numerical integration.
</div>

<div class ="fragment">
<pre>
          QGauss qrule (dim, FIFTH);
        
</pre>
</div>
<div class = "comment">
Tell the finite element object to use our quadrature rule.
</div>

<div class ="fragment">
<pre>
          fe-&gt;attach_quadrature_rule (&qrule);
          
</pre>
</div>
<div class = "comment">
The element Jacobian// quadrature weight at each integration point.   
</div>

<div class ="fragment">
<pre>
          const std::vector&lt;Real&gt;& JxW = fe-&gt;get_JxW();
        
</pre>
</div>
<div class = "comment">
The element shape functions evaluated at the quadrature points.
</div>

<div class ="fragment">
<pre>
          const std::vector&lt;std::vector&lt;Real&gt; &gt;& phi = fe-&gt;get_phi();
          
</pre>
</div>
<div class = "comment">
The element shape function gradients evaluated at the quadrature
points.
</div>

<div class ="fragment">
<pre>
          const std::vector&lt;std::vector&lt;RealGradient&gt; &gt;& dphi = fe-&gt;get_dphi();
        
</pre>
</div>
<div class = "comment">
Now this is slightly different from example 4.
We will not add directly to the overall (PETSc/LASPACK) matrix,
but to the additional matrices "stiffness_mass" and "damping".
The same holds for the right-hand-side vector Fe, which we will
later on store in the additional vector "rhs". 
The zero_matrix is used to explicitly induce the same sparsity
structure in the overall matrix.
see later on. (At least) the mass, and stiffness matrices, however, 
are inherently real.  Therefore, store these as one complex
matrix.  This will definitely save memory.
</div>

<div class ="fragment">
<pre>
          DenseMatrix&lt;Number&gt; Ke, Ce, Me, zero_matrix;
          DenseVector&lt;Number&gt; Fe;
          
</pre>
</div>
<div class = "comment">
This vector will hold the degree of freedom indices for
the element.  These define where in the global system
the element degrees of freedom get mapped.
</div>

<div class ="fragment">
<pre>
          std::vector&lt;unsigned int&gt; dof_indices;
        
</pre>
</div>
<div class = "comment">
Now we will loop over all the elements in the mesh.
We will compute the element matrix and right-hand-side
contribution.


<br><br></div>

<div class ="fragment">
<pre>
          MeshBase::const_element_iterator           el = mesh.local_elements_begin();
          const MeshBase::const_element_iterator end_el = mesh.local_elements_end();
          
          for ( ; el != end_el; ++el)
            {
</pre>
</div>
<div class = "comment">
Start logging the element initialization.
</div>

<div class ="fragment">
<pre>
              START_LOG("elem init","assemble_helmholtz");
              
</pre>
</div>
<div class = "comment">
Store a pointer to the element we are currently
working on.  This allows for nicer syntax later.
</div>

<div class ="fragment">
<pre>
              const Elem* elem = *el;
              
</pre>
</div>
<div class = "comment">
Get the degree of freedom indices for the
current element.  These define where in the global
matrix and right-hand-side this element will
contribute to.
</div>

<div class ="fragment">
<pre>
              dof_map.dof_indices (elem, dof_indices);
        
</pre>
</div>
<div class = "comment">
Compute the element-specific data for the current
element.  This involves computing the location of the
quadrature points (q_point) and the shape functions
(phi, dphi) for the current element.
</div>

<div class ="fragment">
<pre>
              fe-&gt;reinit (elem);
              
</pre>
</div>
<div class = "comment">
Zero & resize the element matrix and right-hand side before
summing them, with different element types in the mesh this
is quite necessary.
</div>

<div class ="fragment">
<pre>
              {
                const unsigned int n_dof_indices = dof_indices.size();
        
        	Ke.resize          (n_dof_indices, n_dof_indices);
        	Ce.resize          (n_dof_indices, n_dof_indices);
        	Me.resize          (n_dof_indices, n_dof_indices);
        	zero_matrix.resize (n_dof_indices, n_dof_indices);
        	Fe.resize          (n_dof_indices);
              }
              
</pre>
</div>
<div class = "comment">
Stop logging the element initialization.