ex7.php

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

</div>

<div class ="fragment">
<pre>
              STOP_LOG("elem init","assemble_helmholtz");
        
</pre>
</div>
<div class = "comment">
Now loop over the quadrature points.  This handles
the numeric integration.
</div>

<div class ="fragment">
<pre>
              START_LOG("stiffness & mass","assemble_helmholtz");
        
              for (unsigned int qp=0; qp&lt;qrule.n_points(); qp++)
        	{
</pre>
</div>
<div class = "comment">
Now we will build the element matrix.  This involves
a double loop to integrate the test funcions (i) against
the trial functions (j).  Note the braces on the rhs
of Ke(i,j): these are quite necessary to finally compute
Real*(Point*Point) = Real, and not something else...
</div>

<div class ="fragment">
<pre>
                  for (unsigned int i=0; i&lt;phi.size(); i++)
        	    for (unsigned int j=0; j&lt;phi.size(); j++)
        	      {
        		Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
        		Me(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
        	      }	  
        	}
        
              STOP_LOG("stiffness & mass","assemble_helmholtz");
        
</pre>
</div>
<div class = "comment">
Now compute the contribution to the element matrix
(due to mixed boundary conditions) if the current
element lies on the boundary. 

<br><br>The following loops over the sides of the element.
If the element has no neighbor on a side then that
side MUST live on a boundary of the domain.
	 

<br><br></div>

<div class ="fragment">
<pre>
              for (unsigned int side=0; side&lt;elem-&gt;n_sides(); side++)
        	if (elem-&gt;neighbor(side) == NULL)
        	  {
        	    START_LOG("damping","assemble_helmholtz");
        	      
</pre>
</div>
<div class = "comment">
Declare a special finite element object for
boundary integration.
</div>

<div class ="fragment">
<pre>
                    AutoPtr&lt;FEBase&gt; fe_face (FEBase::build(dim, fe_type));
        	      
</pre>
</div>
<div class = "comment">
Boundary integration requires one quadraure rule,
with dimensionality one less than the dimensionality
of the element.
</div>

<div class ="fragment">
<pre>
                    QGauss qface(dim-1, SECOND);
        	      
</pre>
</div>
<div class = "comment">
Tell the finte element object to use our
quadrature rule.
</div>

<div class ="fragment">
<pre>
                    fe_face-&gt;attach_quadrature_rule (&qface);
        	      
</pre>
</div>
<div class = "comment">
The value of the shape functions at the quadrature
points.
</div>

<div class ="fragment">
<pre>
                    const std::vector&lt;std::vector&lt;Real&gt; &gt;&  phi_face =
        	      fe_face-&gt;get_phi();
        	      
</pre>
</div>
<div class = "comment">
The Jacobian// Quadrature Weight at the quadrature
points on the face.
</div>

<div class ="fragment">
<pre>
                    const std::vector&lt;Real&gt;& JxW_face = fe_face-&gt;get_JxW();
        	      
</pre>
</div>
<div class = "comment">
Compute the shape function values on the element
face.
</div>

<div class ="fragment">
<pre>
                    fe_face-&gt;reinit(elem, side);
        	      
</pre>
</div>
<div class = "comment">
For the Robin BCs consider a normal admittance an=1
at some parts of the bounfdary
</div>

<div class ="fragment">
<pre>
                    const Real an_value = 1.;
        	      
</pre>
</div>
<div class = "comment">
Loop over the face quadrature points for integration.
</div>

<div class ="fragment">
<pre>
                    for (unsigned int qp=0; qp&lt;qface.n_points(); qp++)
        	      {
        		
</pre>
</div>
<div class = "comment">
Element matrix contributrion due to precribed
admittance boundary conditions.
</div>

<div class ="fragment">
<pre>
                        for (unsigned int i=0; i&lt;phi_face.size(); i++)
        		  for (unsigned int j=0; j&lt;phi_face.size(); j++)
        		    Ce(i,j) += rho*an_value*JxW_face[qp]*phi_face[i][qp]*phi_face[j][qp];
        	      }
        
        	    STOP_LOG("damping","assemble_helmholtz");
        	  }
        
             
</pre>
</div>
<div class = "comment">
Finally, simply add the contributions to the additional
matrices and vector.
</div>

<div class ="fragment">
<pre>
              stiffness.add_matrix      (Ke, dof_indices);
              damping.add_matrix        (Ce, dof_indices);
              mass.add_matrix           (Me, dof_indices);
              
</pre>
</div>
<div class = "comment">
For the overall matrix, explicitly zero the entries where
we added values in the other ones, so that we have 
identical sparsity footprints.
</div>

<div class ="fragment">
<pre>
              matrix.add_matrix(zero_matrix, dof_indices);
        
            } // loop el
        
        
</pre>
</div>
<div class = "comment">
It now remains to compute the rhs. Here, we simply
get the normal velocities values on the boundary from the
mesh data.
</div>

<div class ="fragment">
<pre>
          {
            START_LOG("rhs","assemble_helmholtz");
            
</pre>
</div>
<div class = "comment">
get a reference to the mesh data.
</div>

<div class ="fragment">
<pre>
            const MeshData& mesh_data = es.get_mesh_data();
        
</pre>
</div>
<div class = "comment">
We will now loop over all nodes. In case nodal data 
for a certain node is available in the MeshData, we simply
adopt the corresponding value for the rhs vector.
Note that normal data was given in the mesh data file,
i.e. one value per node
</div>

<div class ="fragment">
<pre>
            assert(mesh_data.n_val_per_node() == 1);
        
            MeshBase::const_node_iterator       node_it  = mesh.nodes_begin();
            const MeshBase::const_node_iterator node_end = mesh.nodes_end();
            
            for ( ; node_it != node_end; ++node_it)
              {
</pre>
</div>
<div class = "comment">
the current node pointer
</div>

<div class ="fragment">
<pre>
                const Node* node = *node_it;
        
</pre>
</div>
<div class = "comment">
check if the mesh data has data for the current node
and do for all components
</div>

<div class ="fragment">
<pre>
                if (mesh_data.has_data(node))
        	  for (unsigned int comp=0; comp&lt;node-&gt;n_comp(0,0); comp++)
        	    {
</pre>
</div>
<div class = "comment">
the dof number
</div>

<div class ="fragment">
<pre>
                      unsigned int dn = node-&gt;dof_number(0,0,comp);
        
</pre>
</div>
<div class = "comment">
set the nodal value
</div>

<div class ="fragment">
<pre>
                      freq_indep_rhs.set(dn, mesh_data.get_data(node)[0]);
        	    }
              }
         
        
            STOP_LOG("rhs","assemble_helmholtz");
          }
        
</pre>
</div>
<div class = "comment">
All done!
</div>

<div class ="fragment">
<pre>
        #endif
        }
        
        
</pre>
</div>
<div class = "comment">
We now define the function which will combine
the previously-assembled mass, stiffness, and
damping matrices into a single system matrix.
</div>

<div class ="fragment">
<pre>
        void add_M_C_K_helmholtz(EquationSystems& es,
        			 const std::string& system_name)
        {
        #ifdef USE_COMPLEX_NUMBERS
        
          START_LOG("init phase","add_M_C_K_helmholtz");
          
</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 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">
Get the frequency, fluid density, and speed of sound
for which we should currently solve
</div>

<div class ="fragment">
<pre>
          const Real frequency = es.parameters.get&lt;Real&gt; ("current frequency");
          const Real rho       = es.parameters.get&lt;Real&gt; ("rho");
          const Real speed     = es.parameters.get&lt;Real&gt; ("wave speed");
          
</pre>
</div>
<div class = "comment">
Compute angular frequency omega and wave number k
</div>

<div class ="fragment">
<pre>
          const Real omega = 2.0*libMesh::pi*frequency;
          const Real k     = omega / speed;
          
</pre>
</div>
<div class = "comment">
Get writable references to the overall matrix and vector, where the 
frequency-dependent system is to be collected
</div>

<div class ="fragment">
<pre>
          SparseMatrix&lt;Number&gt;&  matrix          = *f_system.matrix;
          NumericVector&lt;Number&gt;& rhs             = *f_system.rhs;
          
</pre>
</div>
<div class = "comment">
Get writable references to the frequency-independent matrices
and rhs, though we only need to extract values.  This write access
is necessary, since solver packages have to close the data structure 
before they can extract values for computation.
</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">
form the scaling values for the coming matrix and vector axpy's
</div>

<div class ="fragment">
<pre>
          const Number scale_stiffness (  1., 0.   );
          const Number scale_damping   (  0., omega);
          const Number scale_mass      (-k*k, 0.   );
          const Number scale_rhs       (  0., -(rho*omega));
          
</pre>
</div>
<div class = "comment">
Now simply add the matrices together, store the result
in matrix and rhs.  Clear them first.
</div>

<div class ="fragment">
<pre>
          matrix.close(); matrix.zero ();
          rhs.close();    rhs.zero    ();
          
</pre>
</div>