-#include <base/polynomial_space.h>
-#include <base/parsed_function.h>
-#include <base/smartpointer.h>
-#include <base/convergence_table.h>
+
+
#include <base/quadrature_lib.h>
-#include <base/quadrature_selector.h>
-#include <base/utilities.h>
+#include <base/function.h>
#include <lac/full_matrix.h>
#include <lac/vector.h>
#include <lac/solver_control.h>
-#include <lac/solver_gmres.h>
#include <lac/solver_cg.h>
#include <lac/precondition.h>
-#include <lac/constraint_matrix.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/grid_out.h>
#include <grid/tria_boundary_lib.h>
#include <grid/grid_generator.h>
#include <grid/grid_tools.h>
#include <fstream>
#include <iostream>
-#include <sys/types.h>
-#include <sys/stat.h>
-
-#define deal_II_dimension 3
-
-using std::cout;
-using std::endl;
using namespace dealii;
template <int dim>
class Solution : public Function<dim>
{
-public:
- Solution () : Function<dim>() {}
+ public:
+ Solution () : Function<dim>() {}
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component = 0) const;
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
};
template <int dim>
double Solution<dim>::value (const Point<dim> &p,
- const unsigned int) const
+ const unsigned int) const
{
return sin(numbers::PI * p(0))*cos(numbers::PI * p(1))*exp(p(2));
}
template <int dim>
class RightHandSide : public Function<dim>
{
-public:
- RightHandSide () : Function<dim>() {}
+ public:
+ RightHandSide () : Function<dim>() {}
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
};
template <int dim>
double dPi = numbers::PI;
- // LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu)
+ // LB: u = Delta u - nu D2 u nu - (Grad u nu ) div (nu)
Tensor<2,dim> hessian;
class LaplaceBeltrami
{
public:
- LaplaceBeltrami (Triangulation<dim-1,dim> *tria, Function<dim> &func_data,
- unsigned int fe_degree = 1, unsigned int mapping_degree = 1,
- Function<dim> *pExact = 0);
- // arguments are:
- // triangulation
- // right-hand side
- // fe_degree for solution
- // fe_degree for mapping
- // exact solution is known
+ LaplaceBeltrami ();
+ // arguments are:
+ // triangulation
+ // right-hand side
+ // fe_degree for solution
+ // fe_degree for mapping
+ // exact solution is known
- ~LaplaceBeltrami ();
- void run ();
- double compute_error (VectorTools::NormType norm_type) const;
+ ~LaplaceBeltrami ();
+ void run ();
private:
-
- void setup_system ();
- void assemble_system ();
- void solve ();
- void output_results () const;
+
+ void make_mesh ();
+ void setup_system ();
+ void assemble_system ();
+ void solve ();
+ void output_results () const;
+ void compute_error () const;
- Triangulation<dim-1,dim> *pTria;
- FE_Q<dim-1,dim> fe;
- DoFHandler<dim-1,dim> dh;
- MappingQ<dim-1, dim> mapping;
+ Triangulation<dim-1,dim> triangulation;
+ FE_Q<dim-1,dim> fe;
+ DoFHandler<dim-1,dim> dof_handler;
+ MappingQ<dim-1, dim> mapping;
- ConstraintMatrix matrix_constraints;
+ ConstraintMatrix matrix_constraints;
- SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix;
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
- Vector<double> solution;
- Vector<double> system_rhs;
-
- Function<dim> &rhs_func; // function data
-
- Function<dim> *pExact; // exact solution if provided
-
+ Vector<double> solution;
+ Vector<double> system_rhs;
};
template <int dim>
-LaplaceBeltrami<dim>::LaplaceBeltrami (Triangulation<dim-1,dim> *tria,
- Function<dim> &func_data,
- unsigned int fe_degree,
- unsigned int mapping_degree,
- Function<dim> *pExact
- )
+LaplaceBeltrami<dim>::LaplaceBeltrami ()
:
- fe (fe_degree),
- dh(*tria),
- mapping(mapping_degree),
- rhs_func(func_data)
-{
- pTria = tria;
- this->pExact = pExact;
+ fe (2),
+ dof_handler(triangulation),
+ mapping(2)
+{}
-}
template <int dim>
LaplaceBeltrami<dim>::~LaplaceBeltrami ()
{
- dh.clear ();
+ dof_handler.clear ();
}
+
+template <int dim>
+void LaplaceBeltrami<dim>::make_mesh ()
+{
+ std::map<typename Triangulation<dim-1,dim>::cell_iterator,
+ typename Triangulation<dim,dim>::face_iterator>
+ surface_to_volume_mapping;
+
+ 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);
+
+ static HyperBallBoundary<dim-1,dim> surface_description;
+ triangulation.set_boundary (1, surface_description);
+ triangulation.set_boundary (0, surface_description);
+
+ std::set<unsigned char> boundary_ids;
+ boundary_ids.insert(0);
+
+ GridTools::extract_boundary_mesh (volume_mesh, triangulation,
+ surface_to_volume_mapping,
+ boundary_ids);
+ triangulation.refine_global(3);
+
+ std::cout << "Surface mesh has " << triangulation.n_active_cells()
+ << " cells."
+ << std::endl;
+}
+
+
template <int dim>
void LaplaceBeltrami<dim>::setup_system ()
{
- dh.distribute_dofs (fe);
+ dof_handler.distribute_dofs (fe);
+ std::cout << "Surface mesh has " << dof_handler.n_dofs()
+ << " degrees of freedom."
+ << std::endl;
+
matrix_constraints.clear ();
matrix_constraints.close ();
- CompressedSparsityPattern csp (dh.n_dofs(),
- dh.n_dofs());
+ CompressedSparsityPattern csp (dof_handler.n_dofs(),
+ dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern (dh, csp);
+ DoFTools::make_sparsity_pattern (dof_handler, csp);
matrix_constraints.condense (csp);
sparsity_pattern.copy_from (csp);
system_matrix.reinit (sparsity_pattern);
- solution.reinit (dh.n_dofs());
- system_rhs.reinit (dh.n_dofs());
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
}
template <int dim>
std::vector< double > rhs_values(n_q_points);
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ const RightHandSide<dim> rhs;
+
typename DoFHandler<dim-1,dim>::active_cell_iterator
- cell = dh.begin_active(),
- endc = dh.end();
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
- rhs_func.value_list (fe_values.get_quadrature_points(), rhs_values);
+ rhs.value_list (fe_values.get_quadrature_points(), rhs_values);
for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- for (unsigned int q_point=0; q_point<n_q_points;
- ++q_point)
- {
- cell_matrix(i,j)
- += fe_values.shape_grad(i,q_point) *
- fe_values.shape_grad(j,q_point) *
- fe_values.JxW(q_point);
- }
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ for (unsigned int q_point=0; q_point<n_q_points;
+ ++q_point)
+ {
+ cell_matrix(i,j)
+ += fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point);
+ }
+ }
}
- }
for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
+ {
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- cell_rhs(i) += fe_values.shape_value(i,q_point) *
- rhs_values[q_point]*
- fe_values.JxW(q_point);
- }
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ cell_rhs(i) += fe_values.shape_value(i,q_point) *
+ rhs_values[q_point]*
+ fe_values.JxW(q_point);
+ }
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- system_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- cell_matrix(i,j));
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
- system_rhs(local_dof_indices[i]) += cell_rhs(i);
- }
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
}
- std::map<unsigned int,double> bdy_values;
- VectorTools::interpolate_boundary_values (mapping,dh,
- 0,
- *pExact,
- bdy_values
- );
+ std::map<unsigned int,double> bdy_values;
+ VectorTools::interpolate_boundary_values (mapping,dof_handler,
+ 0,
+ Solution<dim>(),
+ bdy_values
+ );
- MatrixTools::apply_boundary_values (bdy_values,
- system_matrix,
- solution,
- system_rhs,false);
+ MatrixTools::apply_boundary_values (bdy_values,
+ system_matrix,
+ solution,
+ system_rhs,false);
- // condense matrices
+ // condense matrices
matrix_constraints.condense (system_matrix);
matrix_constraints.condense (system_rhs);
}
std::ofstream output (filename.c_str());
DataOut<dim-1,DoFHandler<dim-1,dim> > data_out;
- data_out.attach_dof_handler (dh);
+ data_out.attach_dof_handler (dof_handler);
template <int dim>
-void LaplaceBeltrami<dim>::run ()
-{
-
- setup_system();
- assemble_system();
- solve();
- output_results();
-
+void LaplaceBeltrami<dim>::compute_error () const
+{
+ Vector<float> difference_per_cell (triangulation.n_active_cells());
+ VectorTools::integrate_difference (mapping, dof_handler, solution,
+ Solution<dim>(),
+ difference_per_cell,
+ QGauss<dim-1>(2*fe.degree+1),
+ VectorTools::H1_norm);
+
+
+ std::cout << "H1 error = "
+ << difference_per_cell.l2_norm()
+ << std::endl;
}
-//################################################################################//
-template <int dim>
-double LaplaceBeltrami<dim>::compute_error (VectorTools::NormType norm_type) const
-{
- Assert(pExact != 0, ExcInternalError());
-
- Vector<float> difference_per_cell (pTria->n_active_cells());
- VectorTools::integrate_difference (mapping, dh, solution,
- *pExact, difference_per_cell,
- QGauss<(dim-1)>(2*fe.degree+1),
- norm_type);
- return difference_per_cell.l2_norm();
+template <int dim>
+void LaplaceBeltrami<dim>::run ()
+{
+ make_mesh ();
+ setup_system ();
+ assemble_system ();
+ solve ();
+ output_results ();
+ compute_error ();
}
-
int main ( int argc, char **argv )
{
- std::cout<<std::endl;
- std::cout<<"================================"<<std::endl;
- std::cout<<" LB PROBLEM "<<std::endl;
- std::cout<<"================================"<<std::endl;
- std::cout<<std::endl;
-
-
- Triangulation<deal_II_dimension-1,deal_II_dimension> tria;
-
- // create a mesh consisting on the boundary of the half sphere... thx Seba
- std::map< Triangulation<deal_II_dimension-1,deal_II_dimension>::cell_iterator,
- Triangulation<deal_II_dimension,deal_II_dimension>::face_iterator> surface_to_volume_mapping;
-
- HyperBallBoundary<deal_II_dimension> boundary_description;
- Triangulation<deal_II_dimension> 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);
-
- static HyperBallBoundary<deal_II_dimension-1,deal_II_dimension> surface_description;
- tria.set_boundary (1, surface_description);
- tria.set_boundary (0, surface_description);
-
- std::set<unsigned char> boundary_ids;
- boundary_ids.insert(0);
-
- GridTools::extract_boundary_mesh (volume_mesh, tria,
- surface_to_volume_mapping,
- boundary_ids);
-
-
- tria.refine_global(4);
-
+ try
+ {
+ deallog.depth_console (0);
- RightHandSide<deal_II_dimension> rhs;
- Solution<deal_II_dimension> exact;
-
- LaplaceBeltrami<deal_II_dimension> laplace_beltrami_2d(&tria,rhs,2,2,&exact);
+ 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;
+ }
- laplace_beltrami_2d.run();
- std::cout<<laplace_beltrami_2d.compute_error(VectorTools::H1_norm)<<std::endl;
-
- tria.set_boundary(0);
return 0;
}
-
-
-
-
-
-
-
-
-
-
-