+++ /dev/null
-
-#include "../tests.h"
-#include <deal.II/multigrid/mg_base.h>
-
-/**
- * Generator for system matrix and right hand side.
- */
-class Helmholtz
-{
-public:
- template<int dim, class MATRIX, class VECTOR>
- void build_all(MATRIX &A,
- VECTOR &f,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &quadrature,
- const Function<dim> &rhs);
-
- template<int dim, class MATRIX>
- void build_mgmatrix(MGLevelObject<MATRIX> &A,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &quadrature);
-};
-
-#include "helmholtz1.th"
-#include "helmholtz1mg.th"
-
+++ /dev/null
-
-template<int dim, class MATRIX, class VECTOR>
-void
-Helmholtz::build_all(MATRIX& A,
- VECTOR& f,
- const DoFHandler<dim>& dof,
- const Quadrature<dim>& quadrature,
- const Function<dim>& rhs)
-{
- const unsigned int fe_dofs = dof.get_fe().n_dofs_per_cell();
- FEValues<dim> fe(dof.get_fe(), quadrature,
- UpdateFlags(update_values | update_gradients |
- update_JxW_values | update_q_points));
-
- Vector<double> elvec(fe_dofs);
- FullMatrix<double> elmat(fe_dofs);
- std::vector<unsigned int> indices(fe_dofs);
-
-////////////////////////////////////////////////////
-// Loop for computing the element matrices
-////////////////////////////////////////////////////
-
- for (DoFHandler<dim>::active_cell_iterator c = dof.begin_active()
- ; c != dof.end() ; ++c)
- {
- elmat.clear();
- elvec.clear();
- fe.reinit(c);
-
- // Quadrature loop
-
- for (unsigned k=0;k<quadrature.n_quadrature_points;++k)
- {
- Point<dim> loc = fe.quadrature_point(k);
-
- // Test function loop
- for (unsigned i=0;i<fe_dofs;++i)
- {
- const Point<dim> dv = fe.shape_grad(i,k);
- double v = fe.shape_value(i,k);
-
- elvec(i) += fe.JxW(k) *
- rhs.value(loc) * v;
-
- //Trial function loop
- for (unsigned j=0;j<fe_dofs;++j)
- {
- const Point<dim> du = fe.shape_grad(j,k);
- double u = fe.shape_value(j,k);
-
- elmat(i,j) += fe.JxW(k) *
- (0.1*u*v + du*dv)
- ;
- }
- }
- }
-
- c->get_dof_indices(indices);
- for (unsigned i=0;i<fe_dofs;++i)
- {
- f(indices[i]) += elvec(i);
-
- for (unsigned j=0;j<fe_dofs;++j)
- {
- A.add(indices[i], indices[j], elmat(i,j));
- }
- }
- }
-}
+++ /dev/null
-
-template<int dim, class MATRIX>
-void
-Helmholtz::build_mgmatrix(MGLevelObject<MATRIX>& A,
- const DoFHandler<dim>& dof,
- const Quadrature<dim>& quadrature)
-{
- const unsigned int fe_dofs = dof.get_fe().n_dofs_per_cell();
- FEValues<dim> fe(dof.get_fe(), quadrature,
- UpdateFlags(update_values | update_gradients |
- update_JxW_values | update_q_points));
-
- Vector<double> elvec(fe_dofs);
- FullMatrix<double> elmat(fe_dofs);
- std::vector<unsigned int> indices(fe_dofs);
-
-////////////////////////////////////////////////////
-// Loop for computing the element matrices
-////////////////////////////////////////////////////
-
- DoFHandler<dim>::cell_iterator dc = dof.DoFHandler<dim>::begin();
- for (DoFHandler<dim>::cell_iterator c = dof.begin()
- ; c != dof.end() ; ++c, ++dc)
- {
- elmat.clear();
- elvec.clear();
- fe.reinit(dc);
-
- // Quadrature loop
-
- for (unsigned k=0;k<quadrature.n_quadrature_points;++k)
- {
-// Point<dim> loc = fe.quadrature_point(k);
-
- // Test function loop
- for (unsigned i=0;i<fe_dofs;++i)
- {
- const Point<dim> dv = fe.shape_grad(i,k);
- double v = fe.shape_value(i,k);
-
- //Trial function loop
- for (unsigned j=0;j<fe_dofs;++j)
- {
- const Point<dim> du = fe.shape_grad(j,k);
- double u = fe.shape_value(j,k);
-
- elmat(i,j) += fe.JxW(k) *
- (0.1*u*v + du*dv)
- ;
- }
- }
- }
-
- const unsigned int level = c->level();
-
- c->get_mg_dof_indices(indices);
- for (unsigned i=0;i<fe_dofs;++i)
- {
- for (unsigned j=0;j<fe_dofs;++j)
- {
- A[level].add(indices[i], indices[j], elmat(i,j));
- }
- }
- }
-}