]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
implement first integrators
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 18 Oct 2010 04:50:43 +0000 (04:50 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 18 Oct 2010 04:50:43 +0000 (04:50 +0000)
git-svn-id: https://svn.dealii.org/trunk@22355 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/integrators/laplace.h
deal.II/deal.II/include/integrators/local_integrators.h

index 1863ab571390e7b5b9f02acf8a558b4041ed227c..0b0520932617cc25d9d3dc7b0486a1cc04e40abb 100644 (file)
@@ -33,6 +33,75 @@ namespace LocalIntegrators
  */
   namespace Laplace
   {
+/**
+ * Laplacian in weak form, namely on the cell <i>Z</i> the matrix
+ * \f[
+ * \int_Z \nu \nabla u \cdot \nabla v \, dx.
+ * \f]
+ *
+ * The FiniteElement in <tt>fe</tt> may be scalar or vector valued. In
+ * the latter case, the Laplacian is applied to each component
+ * separately.
+ *
+ * @author Guido Kanschat
+ * @date 2008, 2009, 2010
+ */
+    template<int dim>
+      void cell_matrix (
+       FullMatrix<double>& M,
+       const FEValuesBase<dim>& fe,
+       const double factor = 1.)
+    {
+      const unsigned int n_dofs = fe.dofs_per_cell;
+      const unsigned int n_components = fe.get_fe().n_components();
+      
+      for (unsigned k=0;k<fe.n_quadrature_points;++k)
+       {
+         const double dx = fe.JxW(k) * factor;
+         for (unsigned i=0;i<n_dofs;++i)
+           {
+             for (unsigned j=0;j<n_dofs;++j)
+               for (unsigned int d=0;d<n_components;++d)
+                 M(i,j) += dx *
+                   (fe.shape_grad_component(j,k,d) * fe.shape_grad_component(i,k,d));
+           }
+       }
+    }
+/**
+ * Weak boundary condition of Nitsche type for the Laplacian, namely on the face <i>F<//i> the matrix
+ * @f[
+ * \int_F \Bigl(\gamma u v - \partial_n u v - u \partial_n v\Bigr)\;ds.
+ * @f]
+ *
+ * Here, &gamma; is the <tt>penalty</tt> parameter suitably computed
+ * with compute_penalty().
+ */
+      template <int dim>
+      void nitsche_matrix (
+       FullMatrix<double>& M,
+       const FEValuesBase<dim>& fe,
+       double penalty,
+       double factor = 1.)
+      {
+       const unsigned int n_dofs = fe.dofs_per_cell;
+       const unsigned int n_comp = fe.get_fe().n_components();
+       
+       Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs));
+       Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs));
+       
+       for (unsigned k=0;k<fe.n_quadrature_points;++k)
+         {
+           const double dx = fe.JxW(k) * factor;
+           const Point<dim>& n = fe.normal_vector(k);
+           for (unsigned i=0;i<n_dofs;++i)
+             for (unsigned j=0;j<n_dofs;++j)
+               for (unsigned int d=0;d<n_comp;++d)
+                 M(i,j) += dx *
+                           (fe.shape_value_component(i,k,d) * penalty * fe.shape_value_component(j,k,d)
+                            - (n * fe.shape_grad_component(i,k,d)) * fe.shape_value_component(j,k,d)
+                            - (n * fe.shape_grad_component(j,k,d)) * fe.shape_value_component(i,k,d));
+         }
+      }
   }
 }
 
index da70499a8beeb741637bec0c4e7fdc60dd636af5..2a05d39d40520f6e7a128f36515610b009cdcec0 100644 (file)
@@ -47,7 +47,77 @@ DEAL_II_NAMESPACE_OPEN
  * Functions in this namespace follow a generic signature. In the
  * simplest case, you have two related functions
  * \begin{code}
- * \end{code} 
+ *   template <int dim>
+ *   void
+ *   cell_matrix (
+ *     FullMatrix<double>& M,
+ *     const FEValuesBase<dim>& fe,
+ *     const double factor = 1.);
+ *
+ *   template <int dim>
+ *   void
+ *   cell_residual (
+ *     BlockVector<double>* v,
+ *     const FEValuesBase<dim>& fe,
+ *     const std::vector<Tensor<1,dim> >& input,
+ *     const double factor = 1.);
+ * \end{code}
+ *
+ * There is typically a pair of functions for the same operator, the
+ * function <tt>cell_residual</tt> implementing the mapping of the
+ * operator from the finite element space into its dual, and the
+ * function <tt>cell_matrix</tt> generating the bilinear form
+ * corresponding to the Frechet derivative of <tt>cell_residual</tt>.
+ *
+ * The first argument of these functions is the return type, which is
+ * <ul>
+ * <li> FullMatrix&lt;double&gt; for matrices
+ * <li> BlockVector&ltdouble&gt; for vectors
+ * </li>
+ *
+ * The next argument is the FEValuesBase object representing the
+ * finite element for integration. If the integrated operator maps
+ * from one finite element space into the dual of another (for
+ * instance an off-diagonal matrix in a block system), then first the
+ * FEValuesBase for the trial space and after this the one for the
+ * test space are specified.
+ *
+ * This list is followed by the set of required data in the order
+ * <ol>
+ * <li> Data vectors from finite element functions
+ * <li> Data vectors from other objects
+ * <li> Additional data
+ * <li> A factor which is multiplied with the whole result
+ * </ol>
+ *
+ * <h3>Usage</h3>
+ *
+ * The local integrators can be used wherever a local integration loop
+ * would have been implemented instead. The following example is from
+ * the implementation of a Stokes solver, using
+ * MeshWorker::Assembler::LocalBlocksToGlobalBlocks. The matrices are
+ * <ul>
+ * <li> 0: The vector Laplacian for the velocity (here with a vector valued element)
+ * <li> 1: The divergence matrix
+ * <li> 2: The pressure mass matrix used in the preconditioner
+ * </ul>
+ *
+ * With these matrices, the function called by MeshWorker::loop()
+ * could be written like
+ * \begin{code}
+using namespace ::dealii:: LocalIntegrators;
+
+template <int dim>
+void MatrixIntegrator<dim>::cell(
+  MeshWorker::DoFInfo<dim>& dinfo,
+  typename MeshWorker::IntegrationInfo<dim>& info)
+{
+  Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0));
+  Differential::div_primal_matrix(dinfo.matrix(1,false).matrix, info.fe_values(0), info.fe_values(1));
+  L2::cell_matrix(dinfo.matrix(2,false).matrix, info.fe_values(1));
+}
+ * \end{code}
+ * See step-42 for a worked out example of this code.
  */
 namespace LocalIntegrators
 {

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.