/* to the file deal.II/doc/license.html for the text and */
/* further information on this license. */
+ // @sect3{Include files}
+ // The first few (many?) include
+ // files have already been used in
+ // example 4, so we will
+ // not explain their meaning here
+ // again.
#include <base/quadrature_lib.h>
#include <base/function.h>
-
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_tools.h>
#include <lac/full_matrix.h>
#include <lac/vector.h>
#include <lac/solver_control.h>
#include <lac/precondition.h>
#include <lac/sparse_matrix.h>
#include <lac/compressed_sparsity_pattern.h>
-
-#include <grid/tria.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_boundary_lib.h>
-#include <grid/grid_generator.h>
-#include <grid/grid_tools.h>
-
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
-
#include <fe/fe_values.h>
-
#include <numerics/data_out.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
using namespace dealii;
+ /**
+ @sect3{The <code>LaplaceBeltramiProblem</code> class template}
+
+ This class is extremely similar to the
+ <code>LaplaceProblem</code> class as in
+ example 4.
+ One difference is that now some members
+ will be defined with two template parameters
+ one for the dimension of the mesh @p dim,
+ and the other for the dimension of
+ the embedding space @p spacedim.
+ Now <code>MappingQ</code> appears.
+ */
+template <int spacedim>
+class LaplaceBeltramiProblem
+{
+ private:
+ static const unsigned int dim = spacedim-1;
+
+ public:
+ LaplaceBeltramiProblem (const unsigned degree = 2);
+ void run ();
+
+ private:
+ void make_grid_and_dofs ();
+ void assemble_system ();
+ void solve ();
+ void output_results () const;
+ void compute_error () const;
+
+
+
+ Triangulation<dim,spacedim> triangulation;
+ FE_Q<dim,spacedim> fe;
+ DoFHandler<dim,spacedim> dof_handler;
+ MappingQ<dim, spacedim> mapping;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
template <int dim>
class Solution : public Function<dim>
return return_value;
}
-
+// LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu)
template <int dim>
class RightHandSide : public Function<dim>
{
const unsigned int component = 0) const;
};
-
template <>
double
RightHandSide<3>::value (const Point<3> &p,
{
using numbers::PI;
- // LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu)
+
Tensor<2,3> hessian;
gradient[1] = - PI * sin(PI*p(0))*sin(PI*p(1))*exp(p(2));
gradient[2] = sin(PI*p(0))*cos(PI*p(1))*exp(p(2));
- double curvature;
- Point<3> normal;
- double dLength;
-
- curvature = 3-1; // dim-1
- dLength = sqrt(p(0)*p(0)+p(1)*p(1)+p(2)*p(2));
-
- normal[0] = p(0)/dLength;
- normal[1] = p(1)/dLength;
- normal[2] = p(2)/dLength;
+ Point<3> normal = p;
+ normal /= p.norm();
+
return (-trace(hessian) +
(hessian * normal) * normal +
- (gradient * normal) * curvature);
+ (gradient * normal) * 2.);
}
-
-template <int dim>
-class LaplaceBeltrami
-{
- public:
- LaplaceBeltrami ();
- void run ();
-
- private:
- void make_mesh ();
- void setup_system ();
- void assemble_system ();
- void solve ();
- void output_results () const;
- void compute_error () const;
-
-
-
- Triangulation<dim-1,dim> triangulation;
- FE_Q<dim-1,dim> fe;
- DoFHandler<dim-1,dim> dof_handler;
- MappingQ<dim-1, dim> mapping;
-
- SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix;
-
- Vector<double> solution;
- Vector<double> system_rhs;
-};
-
-
-template <int dim>
-LaplaceBeltrami<dim>::LaplaceBeltrami ()
+template <int spacedim>
+LaplaceBeltramiProblem<spacedim>::LaplaceBeltramiProblem (const unsigned degree)
:
- fe (2),
+ fe (degree),
dof_handler(triangulation),
- mapping(2)
+ mapping(degree)
{}
-template <int dim>
-void LaplaceBeltrami<dim>::make_mesh ()
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::make_grid_and_dofs ()
{
- HyperBallBoundary<dim> boundary_description;
- Triangulation<dim> volume_mesh;
- GridGenerator::half_hyper_ball(volume_mesh);
- volume_mesh.set_boundary (1, boundary_description);
- volume_mesh.set_boundary (0, boundary_description);
- volume_mesh.refine_global (1);
+ Triangulation<spacedim> volume_mesh;
+ GridGenerator::half_hyper_ball(volume_mesh);
- static HyperBallBoundary<dim-1,dim> surface_description;
- triangulation.set_boundary (1, surface_description);
+ static HyperBallBoundary<dim,spacedim> surface_description;
triangulation.set_boundary (0, surface_description);
std::set<unsigned char> boundary_ids;
GridTools::extract_boundary_mesh (volume_mesh, triangulation,
boundary_ids);
- triangulation.refine_global(3);
+ triangulation.refine_global(4);
std::cout << "Surface mesh has " << triangulation.n_active_cells()
<< " cells."
<< std::endl;
-}
-
-template <int dim>
-void LaplaceBeltrami<dim>::setup_system ()
-{
dof_handler.distribute_dofs (fe);
std::cout << "Surface mesh has " << dof_handler.n_dofs()
<< " degrees of freedom."
<< std::endl;
- CompressedSparsityPattern csp (dof_handler.n_dofs(),
- dof_handler.n_dofs());
-
+ CompressedSparsityPattern csp (dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, csp);
-
sparsity_pattern.copy_from (csp);
system_matrix.reinit (sparsity_pattern);
+
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
}
-template <int dim>
-void LaplaceBeltrami<dim>::assemble_system ()
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::assemble_system ()
{
system_matrix = 0;
system_rhs = 0;
- QGauss<dim-1> quadrature_formula(2);
+ QGauss<dim> quadrature_formula(2);
- FEValues<dim-1,dim> fe_values (mapping, fe, quadrature_formula,
- update_values |
- update_gradients |
- update_quadrature_points |
- update_JxW_values);
+ FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula,
+ update_values |
+ update_gradients |
+ update_quadrature_points |
+ update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector< double > rhs_values(n_q_points);
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- const RightHandSide<dim> rhs;
+ const RightHandSide<spacedim> rhs;
- typename DoFHandler<dim-1,dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
-
- for (; cell!=endc; ++cell)
+ for (typename DoFHandler<dim,spacedim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end(); cell!=endc; ++cell)
{
cell_matrix = 0;
cell_rhs = 0;
VectorTools::interpolate_boundary_values (mapping,
dof_handler,
0,
- Solution<dim>(),
+ Solution<spacedim>(),
boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
}
-template <int dim>
-void LaplaceBeltrami<dim>::solve ()
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::solve ()
{
SolverControl solver_control (solution.size(), 1e-7);
SolverCG<> cg (solver_control);
-template <int dim>
-void LaplaceBeltrami<dim>::output_results () const
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::output_results () const
{
- DataOut<dim-1,DoFHandler<dim-1,dim> > data_out;
+ DataOut<dim,DoFHandler<dim,spacedim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution",
- DataOut<dim-1,DoFHandler<dim-1,dim> >::type_dof_data);
+ DataOut<dim,DoFHandler<dim,spacedim> >::type_dof_data);
data_out.build_patches (mapping,
mapping.get_degree());
-template <int dim>
-void LaplaceBeltrami<dim>::compute_error () const
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::compute_error () const
{
Vector<float> difference_per_cell (triangulation.n_active_cells());
VectorTools::integrate_difference (mapping, dof_handler, solution,
- Solution<dim>(),
+ Solution<spacedim>(),
difference_per_cell,
- QGauss<dim-1>(2*fe.degree+1),
+ QGauss<dim>(2*fe.degree+1),
VectorTools::H1_norm);
std::cout << "H1 error = "
-template <int dim>
-void LaplaceBeltrami<dim>::run ()
+template <int spacedim>
+void LaplaceBeltramiProblem<spacedim>::run ()
{
- make_mesh ();
- setup_system ();
+ make_grid_and_dofs();
assemble_system ();
solve ();
output_results ();
int main ( int argc, char **argv )
{
- try
- {
- deallog.depth_console (0);
-
- LaplaceBeltrami<3> laplace_beltrami_2d;
- laplace_beltrami_2d.run();
- }
- catch (std::exception &exc)
- {
- std::cerr << std::endl << std::endl
- << "----------------------------------------------------"
- << std::endl;
- std::cerr << "Exception on processing: " << std::endl
- << exc.what() << std::endl
- << "Aborting!" << std::endl
- << "----------------------------------------------------"
- << std::endl;
-
- return 1;
- }
- catch (...)
- {
- std::cerr << std::endl << std::endl
- << "----------------------------------------------------"
- << std::endl;
- std::cerr << "Unknown exception!" << std::endl
- << "Aborting!" << std::endl
- << "----------------------------------------------------"
- << std::endl;
- return 1;
- }
-
-
+ deallog.depth_console (0);
+ {
+ LaplaceBeltramiProblem<3> laplace_beltrami_2d;
+ laplace_beltrami_2d.run();
+ }
+
return 0;
}
projects / dealii-svn.git / commitdiff